Merge insulin-fix into rc181: a model whose space group gives the same intensities as the written one (same rotations and absences) relabels the data and is refined in its own group

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

image_analysis/structure_refinement/ModelValidation.cpp

@@ -514,6 +514,25 @@ std::string GBytes(size_t bytes) {
 
 }  // namespace
 
+bool SameIntensities(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b) {
+    const gemmi::GroupOps ga = a.operations(), gb = b.operations();
+    const auto rotations = [](const gemmi::GroupOps &g) {
+        std::vector<gemmi::Op::Rot> r;
+        for (const gemmi::Op &op : g.sym_ops)
+            r.push_back(op.rot);
+        std::sort(r.begin(), r.end());
+        return r;
+    };
+    if (rotations(ga) != rotations(gb))
+        return false;
+    for (int h = -6; h <= 6; ++h)
+        for (int k = -6; k <= 6; ++k)
+            for (int l = -6; l <= 6; ++l)
+                if (ga.is_systematically_absent({{h, k, l}}) != gb.is_systematically_absent({{h, k, l}}))
+                    return false;
+    return true;
+}
+
 namespace {
 
 // ValidateAgainstModel() with the rigid body, the structure factors and the maps on the GPU where
@@ -559,17 +578,19 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     }
     result.model_space_group_number = sg->number;
 
-    // If the data was indexed in the enantiomorph of the model's space group (e.g. data P4(1)2(1)2,
-    // model P4(3)2(1)2 - the merged intensities cannot tell them apart), the model's group is a
-    // CANDIDATE for the label the reflections are written under. Only a candidate: this is arithmetic
-    // on two group numbers and says nothing about whether the model belongs to this crystal, and a
-    // model that does not is exactly as capable of rewriting the label as one that does. Whether it is
-    // taken up is settled at the end, on the fit against its own null and on the anomalous map.
+    // If the data were merged in a group the merged intensities cannot tell from the model's - its
+    // enantiomorph (data P4(1)2(1)2, model P4(3)2(1)2), or a group of the same absences (data I23,
+    // model I2(1)3: the search reports the lowest number of such a pair) - the model's group is a
+    // CANDIDATE for the label the reflections are written under, and so for the group the model is
+    // refined in. Only a candidate: this is arithmetic on two groups and says nothing about whether the
+    // model belongs to this crystal, and a model that does not is exactly as capable of rewriting the
+    // label as one that does. Whether it is taken up is settled at the end, on the fit against its own
+    // null and, for an enantiomorph, on the anomalous map.
     //
     // It is tempting to reindex by the change-of-hand operator instead, and that is wrong. The two
-    // groups of an enantiomorphic pair have the same rotation operations - only their translations
-    // differ - so they transform hkl identically, share a reciprocal ASU, and split into Bijvoet
-    // hands identically: the label carries no handedness at all, and nothing about it needs undoing.
+    // groups have the same rotation operations - only their translations differ - so they transform
+    // hkl identically, share a reciprocal ASU, and split into Bijvoet hands identically: the label
+    // carries no handedness at all, and nothing about it needs undoing.
     // What does carry the hand is the indexing the data already have, from the diffraction geometry,
     // and with it the sign of every anomalous difference. The change-of-hand operator is the
     // inversion, so reindexing by it swaps I(+) with I(-) - it does not correct the hand, it flips
@@ -578,20 +599,12 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     // disagreement but manufactures agreement. The anomalous difference map below is the only honest
     // arbiter, and it is used to report the disagreement rather than to bury it.
     const std::vector<MergedReflection> &obs = merged;
-    if (data_space_group && data_space_group->number != sg->number) {
-        if (data_space_group->is_enantiomorphic() && sg->is_enantiomorphic()) {
-            gemmi::GroupOps eops = data_space_group->operations();
-            eops.change_basis_forward(data_space_group->change_of_hand_op());
-            const gemmi::SpaceGroup *enant = gemmi::find_spacegroup_by_ops(eops);
-            if (enant && enant->number == sg->number) {
-                result.model_enantiomorph_candidate = true;
-                logger.Info("Model validation: data space group {} is the enantiomorph of the model {}; "
-                            "the model's group is a candidate for the label, pending the fit - a change "
-                            "of label only, since the two groups index identically and reindexing would "
-                            "flip the anomalous differences",
-                            data_space_group->short_name(), sg->hm);
-            }
-        }
+    if (data_space_group && data_space_group->number != sg->number && SameIntensities(*data_space_group, *sg)) {
+        result.model_enantiomorph_candidate = true;
+        logger.Info("Model validation: data space group {} and the model's {} have the same absences; the "
+                    "model's group is a candidate for the label, pending the fit - a change of label only, "
+                    "since the two groups index identically",
+                    data_space_group->short_name(), sg->hm);
     }
 
     // Resolution limit from the data (the merged set is already resolution-trimmed).
@@ -1774,10 +1787,12 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
     // cannot (measured, inverting the model through the origin moves R-work by less than 1e-4, because
     // |F(h)| of the inverted structure is |F(-h)| = |F(h)| on Friedel-averaged data), so where the
     // anomalous differences do say something they say it alone, and they get a veto.
+    // The veto is the hand's: a pair such as I23 and I2(1)3 is one hand, told apart by the screws in the
+    // model's structure factors, which the fit has already weighed.
     result.adopted_model_enantiomorph = result.model_enantiomorph_candidate && result.model_fits
-                                        && !result.anomalous_hands_disagree;
+                                        && !(sg->is_enantiomorphic() && result.anomalous_hands_disagree);
     if (result.model_enantiomorph_candidate && !result.adopted_model_enantiomorph)
-        logger.Warning("Model validation: the model asserts the enantiomorph {} against the data's {}, "
+        logger.Warning("Model validation: the model asserts its group {} against the data's {}, "
                        "and that assertion is NOT taken up: {}. The reflections are written in the "
                        "group they were merged in",
                        sg->hm, data_space_group != nullptr ? data_space_group->hm : "?",
@@ -2108,7 +2123,7 @@ const gemmi::SpaceGroup *AdoptModelFrame(const ModelValidationResult &validation
     // Bijvoet mates and so change the data - see the note in ValidateAgainstModel.
     if (validation.adopted_model_enantiomorph && validation.model_space_group_number > 0) {
         space_group = gemmi::find_spacegroup_by_number(validation.model_space_group_number);
-        logger.Info("Model validation: the written reflections take the model's enantiomorph, {} ({}), "
+        logger.Info("Model validation: the written reflections take the model's group, {} ({}), "
                     "as a label - no reflection moved",
                     space_group ? space_group->short_name() : "?", validation.model_space_group_number);
     }

image_analysis/structure_refinement/ModelValidation.h

@@ -366,6 +366,11 @@ ModelValidationResult ValidateAgainstModel(const std::vector<MergedReflection> &
                                            const gemmi::Op *twin_law = nullptr,
                                            const gemmi::Structure *placed = nullptr);
 
+// True where merged intensities cannot tell the two groups apart: the same rotations and the same
+// systematic absences. An enantiomorphic pair (P4_1 and P4_3), but also I23 and I2_13 or I222 and
+// I2_12_12_1, where the I-centring already extinguishes what the screw would.
+bool SameIntensities(const gemmi::SpaceGroup &a, const gemmi::SpaceGroup &b);
+
 // Reindex `merged` into the frame ValidateAgainstModel reported, so the reflection files that are
 // written describe the same indexing as the R-factors and the maps. Returns the space group they are
 // now in - the model's, where its enantiomorph was adopted, and `data_space_group` otherwise

rugnux/ResultReport.cpp

@@ -2878,9 +2878,9 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
                              " says how - so the reflection files, the R-factors and the maps all describe one\n"
                              " indexing."));
             if (mv.adopted_model_enantiomorph)
-                Add(s, Prose(" The written reflections carry the model's enantiomorph as their space group. That is\n"
-                             " a change of label and nothing else: the two groups of an enantiomorphic pair have\n"
-                             " the same rotation operations, so no reflection moved."));
+                Add(s, Prose(" The written reflections carry the model's space group - its enantiomorph, or a group\n"
+                             " of the same absences - as their own. That is a change of label and nothing else: the\n"
+                             " two groups have the same rotation operations, so no reflection moved."));
             else if (mv.model_enantiomorph_candidate)
                 Add(s, Prose(" The model asserts the other enantiomorph of the written space group, and that\n"
                              " assertion was NOT taken up: merged intensities cannot check it, and the anomalous\n"

rugnux/RugnuxRefine.cpp

@@ -8,6 +8,7 @@
 #include "../common/UnitCell.h"
 #include "../image_analysis/scale_merge/FrenchWilson.h"
 #include "../image_analysis/scale_merge/RfreeFlags.h"
+#include "../image_analysis/structure_refinement/ModelValidation.h"   // SameIntensities
 
 RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &model, double d_min_limit,
                                const RefinementSettings &settings, Logger &logger) {
@@ -19,6 +20,14 @@ RefineOnlyResult RunRefineOnly(RugnuxReflections data, const gemmi::Structure &m
         out.refinement.failure_reason = "unknown space group " + data.space_group;
         return out;
     }
+    // A file labelled with the other group of a pair the intensities cannot separate (I23 for an I2(1)3
+    // model, P4(1) for a P4(3) one) describes the model's group just as well, and the screws are the model's.
+    if (const gemmi::SpaceGroup *msg = model.find_spacegroup();
+        msg != nullptr && msg->number != sg->number && SameIntensities(*sg, *msg)) {
+        logger.Info("The model names space group {}, the data {}, which the intensities cannot tell apart: "
+                    "refined in the model's", msg->xhm(), sg->xhm());
+        sg = msg;
+    }
     out.space_group = sg->xhm();
     out.wavelength_A = data.wavelength_A > 0 ? data.wavelength_A : 1.0;
     std::vector<MergedReflection> &refl = data.reflections;

rugnux/RugnuxScaleMerge.cpp

@@ -5101,17 +5101,37 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
             // A cut (the broker's deadline or cancel, Ctrl+C) cuts a running refinement short, keeping what it has
             settings.cancel = cancelled_.RunFlag();
             const PlacedModelOutput &p = *placed_model_out;
+            // In the group the model was validated in: its own, where that has the rotations of the group
+            // written - the screws are then the model's, as its fit already took them, with the reflections
+            // that group extinguishes left out. Refined in the written group instead, a model whose screws
+            // the merge did not decide (I2(1)3 written I23, P6_3 written P6 with no 00l measured), or whose
+            // screws lie on other axes (P2_12_12 written, the model's P2_12_2_1 on the same cell), packs its
+            // copies wrongly and refines from an R-free near 0.5.
+            gemmi::SpaceGroup refine_sg = p.space_group;
+            std::vector<MergedReflection> refine_merged;
+            if (const gemmi::SpaceGroup *own = p.model->find_spacegroup();
+                own != nullptr && own->xhm() != refine_sg.xhm() && own->centring_type() == refine_sg.centring_type()
+                && SortedRotations(own->operations()) == SortedRotations(refine_sg.operations())) {
+                logger.Info("Refinement: in the model's {}, the group it was validated in, not the written {}",
+                            own->xhm(), refine_sg.xhm());
+                refine_sg = *own;
+                const gemmi::GroupOps ops = own->operations();
+                for (const MergedReflection &r : *placed_merged)
+                    if (!ops.is_systematically_absent({{r.h, r.k, r.l}}))
+                        refine_merged.push_back(r);
+            }
             // processing-inert: a refinement that fails is reported, and the run goes on to its other outputs
             RefinementResult refined;
             try {
-                refined = RefineModel(*p.model, *placed_merged, gemmi::UnitCell(p.cell), p.space_group,
-                                      experiment_.GetWavelength_A(), logger, settings);
+                refined = RefineModel(*p.model, refine_merged.empty() ? *placed_merged : refine_merged,
+                                      gemmi::UnitCell(p.cell), refine_sg, experiment_.GetWavelength_A(), logger,
+                                      settings);
             } catch (const std::exception &e) {
                 refined = RefinementResult();
                 refined.failure_reason = std::string("refinement failed: ") + e.what();
             }
             if (refined.ok) {
-                result.outputs->refined_model = PlacedModelOutput{refined.model, p.cell, p.space_group};
+                result.outputs->refined_model = PlacedModelOutput{refined.model, p.cell, refine_sg};
                 result.outputs->refined_maps = refined.maps;
                 result.outputs->blobs = refined.blobs;
                 result.outputs->solvent = refined.solvent;

tests/ModelValidationTest.cpp

@@ -972,6 +972,22 @@ TEST_CASE("ModelValidation_WrongModelInAGroupOffTheDataMetric", "[ModelValidatio
     std::filesystem::remove(wrong_path);
 }
 
+// The groups a model may relabel the data into: those the merged intensities cannot separate from the
+// data's - an enantiomorph, or a group whose screws the centring already hides - and no other.
+TEST_CASE("ModelValidation_SameIntensities", "[ModelValidation]") {
+    const auto same = [](const char *a, const char *b) {
+        return SameIntensities(*gemmi::find_spacegroup_by_name(a), *gemmi::find_spacegroup_by_name(b));
+    };
+    CHECK(same("I 2 3", "I 21 3"));
+    CHECK(same("I 2 2 2", "I 21 21 21"));
+    CHECK(same("P 41 21 2", "P 43 21 2"));
+    CHECK(same("P 61", "P 65"));
+    CHECK(!same("P 2 2 2", "P 21 21 21"));
+    CHECK(!same("P 61", "P 62"));
+    CHECK(!same("P 2 3", "I 2 3"));
+    CHECK(!same("I 4 3 2", "I 2 3"));
+}
+
 // A model of something else entirely, scaled to the data, shares their falloff across every shell, and
 // by Fisher's transform on independent reflections that alone reads as a significant correlation. Against
 // the null of the data at h with the model at h + d - which shares the falloff too - it must not.