diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index 50afeb0f2..c2b2ca749 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/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 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 &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 &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 &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 &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); } diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index 5bb6d3043..a42be37d3 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -366,6 +366,11 @@ ModelValidationResult ValidateAgainstModel(const std::vector & 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 diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 179004b29..424856e7b 100644 --- a/rugnux/ResultReport.cpp +++ b/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" diff --git a/rugnux/RugnuxRefine.cpp b/rugnux/RugnuxRefine.cpp index d232a8074..fd8e51d90 100644 --- a/rugnux/RugnuxRefine.cpp +++ b/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 &refl = data.reflections; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 4bacea24e..978469bb3 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/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 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; diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index d0d71b1dd..7b622a91d 100644 --- a/tests/ModelValidationTest.cpp +++ b/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.