From d9b025b5af1c0e093e80be37727eda88e7ef002b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 12:21:42 +0100 Subject: [PATCH 1/3] RESULT-CHANGING: symmetry hypotheses of a refused promotion carried to the model and decided there Where the space-group search refuses a higher point group G over the adopted A, every other index-2 subgroup S of G perfectly twinned by a law of its coset merges into G exactly as A does, so no intensity statistic separates them; the model can. Following Zwart et al. (2008) and Lebedev & Isupov (2014): refine in the candidate groups and compare. - SiblingTwinSubgroups (TwinningAnalysis): the siblings of A in G with their coset laws; Catch2 test. - Report: SYMMETRY_HYPOTHESES / SYMMETRY_HYPOTHESIS_n (adopted, adopted twinned, each sibling twinned, and whether the model tested it) on every run with a refusal. Report only without a model. - --model: where the model's own group has a sibling's rotations on the data's axes, that sibling's merge is the twin hypothesis the model is scored on (twin target, twin-closed free set from the lattice-holohedry orbits), and the placement is refined in the model's own group against the model-choice file (absences of that group dropped). - Best two into refinement: where the model fits the answer's merge and the twin hypothesis' within 0.03 in R-free, both are refined (each with its own null-law twin decision) and the rival replaces the validation's choice only where it refines better by 0.01. MODEL_CHOICE_DECIDED_BY, MODEL_CHOICE_TWIN_LAW, RIVAL_REFINE_*. - MR: sibling merges are searched beside the answer's (every sibling for a model without a group of its own) and chosen where the twinned R-free beats the answer's best by R_FREE_MARGIN. A model choice MR makes (another group or a twin hypothesis) is now refined in that group against its file. Gated on the existing refusal: no new thresholds; sets without a refused promotion run as before. On a perfectly merohedrally twinned P6x crystal with three NCS copies and strong tNCS (an NCS 2-fold parallel to the twin axis), adopted as P3x21: --model now places and refines in P6x twinned by the in-plane law, where the refinement in P3x21 diverged. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/ACKNOWLEDGEMENT.md | 11 + docs/RUGNUX_REPORT.md | 23 ++ .../scale_merge/TwinningAnalysis.cpp | 19 ++ image_analysis/scale_merge/TwinningAnalysis.h | 8 + .../MolecularReplacement.cpp | 45 ++- .../MolecularReplacement.h | 3 + rugnux/ResultReport.cpp | 42 +++ rugnux/Rugnux.h | 29 ++ rugnux/RugnuxScaleMerge.cpp | 318 +++++++++++++++--- tests/TwinningAnalysisTest.cpp | 34 ++ 10 files changed, 492 insertions(+), 40 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index ddbb28396..ede0116af 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -229,6 +229,17 @@ H = |I(h) − I(Th)|/(I(h) + I(Th)) of the pairs a law relates is uniform on [0, gives the fraction. T. O. Yeates, "Detecting and overcoming crystal twinning" (1997), Methods Enzymol. 276, 344-358 [doi:10.1016/S0076-6879(97)76068-3](https://doi.org/10.1016/S0076-6879%2897%2976068-3). +**The symmetry hypotheses a refused promotion leaves** - the adopted group and every other index-2 subgroup +of the refused group perfectly twinned by a law of its coset, all merging alike, carried to the model and +decided by placing and refining it in each - follow the advice to refine in the candidate groups and compare, +rather than decide from intensity statistics, that pseudo-symmetry and twinning call for: +P. H. Zwart, R. W. Grosse-Kunstleve, A. A. Lebedev, G. N. Murshudov and P. D. Adams, "Surprises and pitfalls +arising from (pseudo)symmetry" (2008), Acta Cryst. D64, 99-107 +[doi:10.1107/S090744490705531X](https://doi.org/10.1107/S090744490705531X); A. A. Lebedev and M. N. Isupov, +"Space-group and origin ambiguity in macromolecular structures with pseudo-symmetry and its treatment with the +program Zanuda" (2014), Acta Cryst. D70, 2430-2443 +[doi:10.1107/S1399004714014795](https://doi.org/10.1107/S1399004714014795). + **Twin refinement** (§15.12) refines the model against Yeates's twinned intensity, (1−α) I(h) + α I(Th), with the twin fraction refined beside it as SHELXL refines BASF: R. Herbst-Irmer and G. M. Sheldrick, "Refinement of twinned structures with SHELXL97" (1998), Acta Cryst. B54, 443-449 diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 3bbbc2928..4820546b9 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -996,6 +996,29 @@ twin hypothesis' merge is also made where a metric re-ask promoted the indexed l the model's group has fewer rotations than the adopted one; made for that reason alone, it is written and reported only where the model chooses it. +**Symmetry hypotheses of a refused promotion.** Where the search refused a higher point group, the +adopted group is one of several hypotheses the merged intensities cannot separate: the adopted group +untwinned or twinned by a law of the refused group outside it, and each *other* index-2 subgroup of the +refused group perfectly twinned by a law of its coset - each of them merges into the refused group alike. +The report lists them, on every run, as `SYMMETRY_HYPOTHESES=` (their number) and +`SYMMETRY_HYPOTHESIS_= | | ADOPTED / ADOPTED_TWINNED / SIBLING_TWINNED | +TESTED_BY_MODEL / UNTESTED`; a sibling's group is named by its rotations alone (symmorphic), its screws +left to the model. They cost nothing without a model. With `--model`, where the model's own group is a +sibling's (it has the sibling's rotations on the data's axes), that sibling's merge is the twin hypothesis +the model is scored on, exactly as a subgroup's above, and the model is refined in its own group there. +With a model that has no group of its own (a prediction, or `--mr` with coordinates without a cell), every +sibling's merge is searched by molecular replacement beside the answer's. + +Where the model fits both the answer's merge and the twin hypothesis' within 0.03 in R-free before +refinement, **both are refined**, each with its own twin decision (below), and the refined R-free decides +which the model chooses: the rival replaces the validation's choice only where it refines better by 0.01, +what two refinements of one model differ by from their paths alone. `MODEL_CHOICE_DECIDED_BY=` is +`VALIDATION` or `REFINEMENT`; the rival's `RIVAL_REFINE_SPACE_GROUP`, `RIVAL_REFINE_MERGE`, +`RIVAL_REFINE_R_WORK`, `RIVAL_REFINE_R_FREE`, `RIVAL_REFINE_TWIN` (with `RIVAL_REFINE_TWIN_LAW`) and +`RIVAL_REFINE_REPLACED` stand beside the `REFINE_` keys of the one kept. `MODEL_CHOICE_TWIN_LAW=` is the law +the model's choice file carries (`NONE` without one). The model's choice is refined against its own files +by either route, molecular replacement's included. + The null's own numbers — `MODEL_FIT_STATISTIC`, `MODEL_FIT_VALUE`, the `MODEL_FIT_NULL_*` keys and the `MODEL_INDEXING_MARGIN*` keys — are written with `--developer`; the default report carries the verdict (`MODEL_FIT`, `MODEL_FIT_SIGMA`, `MODEL_DECISIONS_TAKEN` and the two decisions). diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index a5cfae99a..79af2e7fb 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -60,6 +60,25 @@ namespace { } } +std::vector SiblingTwinSubgroups(const gemmi::SpaceGroup& adopted, const gemmi::SpaceGroup& refused) { + const auto rotations = [](const gemmi::SpaceGroup& g) { + std::vector r; + for (const auto& op : g.operations().sym_ops) + r.push_back(op.rot); + std::sort(r.begin(), r.end()); + return r; + }; + const auto a = rotations(adopted); + const auto g = rotations(refused); + if (2 * a.size() != g.size() || !std::includes(g.begin(), g.end(), a.begin(), a.end())) + return {}; + std::vector out; + for (auto& t : IndexTwoTwinSubgroups(refused)) + if (rotations(t.subgroup) != a) + out.push_back(std::move(t)); + return out; +} + bool TwinLawPossible(const gemmi::UnitCell& cell, const gemmi::SpaceGroup& space_group, double max_obliquity_deg) { return !gemmi::find_twin_laws(cell, &space_group, max_obliquity_deg, /*all_ops=*/false).empty(); } diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index 7a2d74835..c4d40b2e0 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -126,6 +126,14 @@ struct TwinSubgroup { }; std::vector IndexTwoTwinSubgroups(const gemmi::SpaceGroup& group); +// Following Zwart et al. (2008) Acta Cryst. D64, 99-107, and Lebedev & Isupov (2014) Acta Cryst. D70, 2430-2443: +// the other perfect-twin hypotheses a refused promotion leaves. Where `adopted` is an index-2 subgroup of `refused`, every OTHER index-2 subgroup of `refused` +// (IndexTwoTwinSubgroups) with the two-folds of its coset as laws. Each, perfectly twinned by its law, +// merges into `refused` exactly as `adopted` does, so the merged intensities hold nothing that tells +// them apart - only a model does. Both groups in the same basis; the subgroups keep `refused`'s +// translations. Empty where `adopted` is not an index-2 subgroup of `refused`. +std::vector SiblingTwinSubgroups(const gemmi::SpaceGroup& adopted, const gemmi::SpaceGroup& refused); + // The twin fraction a twin law implies on a merge in a group the law is not part of, from the // disagreement of the reflection pairs it relates: H = |I(h) - I(Th)| / (I(h) + I(Th)) over the // acentric pairs, uniform on [0,1] for untwinned Wilson intensities and on [0, 1-2a] for a twin of diff --git a/image_analysis/structure_refinement/MolecularReplacement.cpp b/image_analysis/structure_refinement/MolecularReplacement.cpp index 816f0dc18..dc18038e8 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.cpp +++ b/image_analysis/structure_refinement/MolecularReplacement.cpp @@ -3482,7 +3482,48 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec // merge cannot hold the model in any of its groups. search_merge(0); const int answer = choose_in(0); - if (answer >= 0) { + // A sibling hypothesis' merge (MRMerge::sibling): the same merged intensities as the answer's could come from + // either, so it is searched whatever the answer's merge holds, and chosen where its own placement fits its data + // as twinned by its law better by R_FREE_MARGIN in R-free than the answer's choice fits the answer's merge - or + // where the answer's merge chose nothing - and passes the same choice as any merge's. + const bool siblings = merges.size() > 1 && merges[1].sibling; + if (siblings) { + result.chosen = answer >= 0 ? answer : -1; + const std::string decision = result.decision; + double answer_r_free = INFINITY; // the answer's merge's best placement, chosen or not + for (const MRCandidate &o : result.candidates) + if (o.merge == 0 && o.verdict == "PLACED" && o.r_free < answer_r_free) + answer_r_free = o.r_free; + int best = -1; + bool any_placed = false; + for (size_t m = 1; m < merges.size(); m++) { + if (!merges[m].sibling || !merges[m].twin_law) + continue; + search_merge(m); + const int h = choose_in(static_cast(m)); + any_placed = any_placed || h != -1; + if (h < 0) + continue; + MRCandidate &c = result.candidates[h]; + const PlacementScore t = score_placement(*c.placed, c.space_group, *merges[m].merged, result.d_min_score, + nthreads, &scorers, nullptr, &*merges[m].twin_law); + c.r_work_twin = t.r_work_twin; + c.r_free_twin = t.r_free_twin; + c.twin_fraction = t.twin_fraction; + logger.Info("Molecular replacement: sibling hypothesis {} twinned by {}: {} copies, R-free {:.3f} at a " + "fraction {:.2f}, against {:.3f} in the answer's merge", + c.space_group.xhm(), merges[m].twin_law->as_hkl().triplet('h'), c.copies, t.r_free_twin, + t.twin_fraction, answer_r_free); + if (best < 0 || t.r_free_twin < result.candidates[best].r_free_twin) + best = h; + } + if (best >= 0 && result.candidates[best].r_free_twin <= answer_r_free - R_FREE_MARGIN) + result.chosen = best; + if (result.chosen >= 0 || answer == -2) + result.decision = decision; + else if (answer == -1 && !any_placed) + result.decision = "NOT_PLACED: no group of any merge holds a placed model"; + } else if (answer >= 0) { result.chosen = answer; } else if (answer == -1 && merges.size() > 1 && cannot_hold(0)) { search_merge(1); @@ -3504,7 +3545,7 @@ MolecularReplacementResult search_model(const SearchModel &model, const std::vec // copies, nor the fraction, was fitted to - and passes the same choice as any merge's (its null, and a // lead over the subgroup's other groups). if (result.chosen >= 0 && result.candidates[result.chosen].merge == 0 && merges.size() > 1 - && merges[1].twin_law) { + && merges[1].twin_law && !siblings) { const MRCandidate g = result.candidates[result.chosen]; const MRMerge &tm = merges[1]; const gemmi::UnitCell h_cell(tm.cell); diff --git a/image_analysis/structure_refinement/MolecularReplacement.h b/image_analysis/structure_refinement/MolecularReplacement.h index d56d402f0..3784b07d9 100644 --- a/image_analysis/structure_refinement/MolecularReplacement.h +++ b/image_analysis/structure_refinement/MolecularReplacement.h @@ -310,6 +310,9 @@ struct MRMerge { UnitCell cell; gemmi::SpaceGroup space_group; // the group it was merged in std::optional twin_law; // the twin hypothesis' merge: the law it carries + // The twin hypothesis is not a subgroup of the answer's group but another index-2 subgroup of a group the + // space-group search refused (SiblingTwinSubgroups): it merges as the answer does, so it is searched beside it. + bool sibling = false; }; // One space group's best placement: its first copy and every further one the search added. diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 179004b29..ba9912d4b 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1189,6 +1189,25 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyText("SPACE_GROUP_REFUSED_POINT_GROUP", result.space_group_search->refused_point_group_hm)); Add(s, KeyText("SPACE_GROUP_REFUSED_REASON", result.space_group_search->refused_reason)); + // The hypotheses the refusal leaves, which no intensity statistic separates (each, perfectly + // twinned by its law, merges into the refused group alike) - for a model to decide. + if (!result.symmetry_hypotheses.empty()) { + Add(s, KeyInt("SYMMETRY_HYPOTHESES", static_cast(result.symmetry_hypotheses.size()))); + for (size_t i = 0; i < result.symmetry_hypotheses.size(); i++) { + const auto &h = result.symmetry_hypotheses[i]; + Add(s, KeyText(fmt::format("SYMMETRY_HYPOTHESIS_{}", i + 1).c_str(), + fmt::format("{} | {} | {} | {}", h.space_group, + h.twin_law.empty() ? "NONE" : h.twin_law, h.kind, + h.tested_by_model ? "TESTED_BY_MODEL" : "UNTESTED"))); + } + Add(s, Prose(fmt::format( + " The refused {} leaves these symmetry hypotheses: the adopted group, untwinned and twinned\n" + " by a law of the refused group outside it, and each other index-2 subgroup of the refused\n" + " group perfectly twinned by a law of its coset. All of them merge into the refused group\n" + " alike, so the intensities cannot choose; a model can (--model, --mr). A sibling's group\n" + " is named without screws, which the model or the axial rows decide.", + result.space_group_search->refused_point_group_hm))); + } } else { Add(s, KeyText("SPACE_GROUP_REFUSED_POINT_GROUP", "NONE", true)); } @@ -2743,6 +2762,29 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, if (!result.model_prior_file_suffix.empty()) { Add(s, KeyText("MODEL_CHOICE_FILE_PREFIX", output_prefix + result.model_prior_file_suffix)); Add(s, KeyText("MODEL_CHOICE_SPACE_GROUP", result.model_prior_file_space_group)); + Add(s, KeyText("MODEL_CHOICE_TWIN_LAW", result.model_prior_file_twin_law.empty() + ? std::string("NONE") : result.model_prior_file_twin_law)); + } + if (!result.model_choice_decided_by.empty()) + Add(s, KeyEnum("MODEL_CHOICE_DECIDED_BY", result.model_choice_decided_by, {"VALIDATION", "REFINEMENT"})); + // Two symmetry hypotheses the model fitted alike, both refined: the rival's numbers, beside the + // REFINE_ keys of the one kept. + if (result.rival_refinement) { + const auto &rr = *result.rival_refinement; + Add(s, KeyText("RIVAL_REFINE_SPACE_GROUP", rr.space_group)); + Add(s, KeyEnum("RIVAL_REFINE_MERGE", rr.merge, {"MERGE", "TWIN_SUBGROUP_MERGE"})); + Add(s, KeyReal("RIVAL_REFINE_R_WORK", rr.r_work, "{:.4f}")); + Add(s, KeyReal("RIVAL_REFINE_R_FREE", rr.r_free, "{:.4f}")); + Add(s, KeyText("RIVAL_REFINE_TWIN", rr.twin.empty() ? std::string("NONE") : rr.twin)); + if (!rr.twin_law.empty()) + Add(s, KeyText("RIVAL_REFINE_TWIN_LAW", rr.twin_law)); + Add(s, KeyBool("RIVAL_REFINE_REPLACED", rr.replaced)); + Add(s, Prose(fmt::format( + " The model fitted the answer's merge and the twin hypothesis' within {:.2f} in R-free\n" + " before refinement, so both were refined. {} refined to R-free {:.4f}, the rival {} to\n" + " {:.4f}; the rival {} the model's choice (it has to refine better by {:.2f}).", + 0.03, rr.chosen_space_group, rr.chosen_r_free, rr.space_group, rr.r_free, + rr.replaced ? "became" : "did not replace", 0.01))); } // The twin hypothesis (section 4): the model scored against that merge too, and which of // the two the maps and the placed model were made from. The keys above are the answer's. diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2db093e7c..08f460ef7 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -751,6 +751,35 @@ struct ProcessResult { // (_mr), and that group; empty where nothing was. std::string model_prior_file_suffix; std::string model_prior_file_space_group; + std::string model_prior_file_twin_law; // the law that file carries; empty where it carries none + // VALIDATION where the model's choice was the unrefined fit's alone, REFINEMENT where two hypotheses were + // refined and the refined R-free decided (rival_refinement). Empty without a model. + std::string model_choice_decided_by; + // The other symmetry hypothesis refined beside the model's choice (RivalPlacement in RugnuxScaleMerge.cpp): + // where the model fitted the answer's merge and the twin hypothesis' within a margin, both were refined. + struct RivalRefinement { + std::string space_group; // the rival's, as refined + std::string merge; // MERGE or TWIN_SUBGROUP_MERGE + double r_work = NAN, r_free = NAN; + std::string twin, twin_law; // its refinement's twin decision + std::string chosen_space_group; // the validation's choice, as refined + double chosen_r_free = NAN; + bool replaced = false; // the rival refined better and became the model's choice + }; + std::optional rival_refinement; + + // The symmetry hypotheses the merged intensities cannot separate, where the search refused a higher + // point group: the adopted group, and each other index-2 subgroup of the refused group perfectly + // twinned by a law of its coset (SiblingTwinSubgroups) - each merges into the refused group alike, so + // only a model decides. The groups are symmorphic where the screws are open; in the written axes. + // Empty where nothing was refused. + struct SymmetryHypothesis { + std::string space_group; + std::string twin_law; // empty: untwinned + std::string kind; // ADOPTED, ADOPTED_TWINNED, SIBLING_TWINNED + bool tested_by_model = false; // the model was scored on this hypothesis' merge + }; + std::vector symmetry_hypotheses; // Space group used for the final re-scale/merge and written to the master file: determined by the // search when the user did not fix one, otherwise the fixed group. The whole group, so it keeps diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 4bacea24e..78eec1738 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -161,6 +162,17 @@ namespace { // of this and what the geometry asks of the centre (BeamCenterNeed_pxl). constexpr float MAX_BEAM_CENTER_SIGMA_PXL = 1.0f; + // Two symmetry hypotheses a model fits this close in R-free, unrefined - the answer's merge and the twin + // hypothesis' - are both refined, and decided on the refined R-free: molecular replacement's margin for + // telling two placements apart (R_FREE_MARGIN, MolecularReplacement.cpp). Following Zwart et al. (2008) + // Acta Cryst. D64, 99-107, and Lebedev & Isupov (2014) Acta Cryst. D70, 2430-2443: refine in the + // candidate groups and compare. + constexpr double HYPOTHESIS_R_FREE_MARGIN = 0.03; + // What two refinements of one model on one merge differ by in R-free from their paths alone (a twinned + // and an untwinned arm of the same model that the null puts level, ~0.01): the rival refined has to + // beat the hypothesis the validation chose by more than this to replace it. + constexpr double REFINEMENT_PATH_NOISE = 0.01; + // The twin-immune zone of the operators `higher` adds over `lower`, read on a P1 merge, as a // verdict: +1 the added operators are real (the zone is centric with margin), -1 they are a twin // law or a pseudo-symmetry (acentric with margin), 0 undecided - fewer than two orders apart @@ -327,6 +339,10 @@ namespace { gemmi::SpaceGroup subgroup; gemmi::Op law; std::vector others; + // Not a subgroup of the adopted group but another index-2 subgroup of the group the search refused + // (SiblingTwinSubgroups): perfectly twinned by `law`, it merges into the refused group as the adopted + // one does. + bool sibling = false; }; std::optional ChooseTwinSubgroup(const gemmi::SpaceGroup &group, @@ -376,6 +392,36 @@ namespace { return best; } + // The perfect-twin hypotheses beside the adopted group that a refused promotion leaves: each other + // index-2 subgroup of the refused group (SiblingTwinSubgroups), with the law of its coset that agrees + // worst with the intensities (lowest Stage A CC), as ChooseTwinSubgroup takes it. Given the rotations of + // a model's group, only the sibling with exactly those. Empty where nothing was refused. + std::vector SiblingTwinChoices(const gemmi::SpaceGroup &adopted, + const SearchSpaceGroupResult &search, + const std::vector> &model = {}) { + std::vector out; + if (!search.refused_point_group_representative) + return out; + const auto law_cc = [&](const gemmi::Op &law) { + const std::string t = law.as_hkl().triplet('h'); + for (const auto &s : search.operator_scores) + if (s.op_triplet_hkl == t) + return s.cc; + return std::numeric_limits::infinity(); + }; + for (const auto &t : SiblingTwinSubgroups(adopted, *search.refused_point_group_representative)) { + if (!model.empty() + && std::find(model.begin(), model.end(), SortedRotations(t.subgroup.operations())) == model.end()) + continue; + const gemmi::Op *law = &t.laws.front(); + for (const auto &l : t.laws) + if (law_cc(l) < law_cc(*law)) + law = &l; + out.push_back(TwinSubgroupChoice{t.subgroup, *law, {}, true}); + } + return out; + } + // A twin-law choice on other axes: the group by SpaceGroupInBasis, the law conjugated by the // Miller-index operator the change of basis moves the reflections by. False where no setting in // GEMMI's table names the moved group. @@ -1987,6 +2033,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // The model's own subgroup, where it is not the one above (which the data chose): a merge of its own, // made only for the model to test. std::optional model_twin_choice; + // The perfect-twin hypotheses a refused promotion leaves beside the adopted group, in the scaler's frame + // (SiblingTwinChoices): reported, and the model's own among them tested as model_twin_choice. + std::vector sibling_choices; + std::string sibling_for_model; // the one of them whose merge was made for the model, by name (scaler's frame) // The --model's kind, from the file alone: a prediction has no placement of its own and goes straight to // molecular replacement; an experimental structure's group, where it has fewer rotations than the merge's, // asks for the twin hypothesis' merge. @@ -3572,9 +3622,22 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // A model can test the twin hypothesis itself (the twin validation, molecular replacement's twin // comparison), so with one the merge in the subgroup is made wherever a metric re-ask promoted from the // indexed lattice's group, or the model is in a group with fewer rotations than the adopted one. + // The perfect-twin hypotheses a refused promotion leaves beside the adopted group (SiblingTwinChoices), + // reported on every run; and the one among them a model's own group is in - the model tests it as it + // tests a subgroup's. + if (!geometry_prepass && result.space_group && result.consensus_cell) + sibling_choices = SiblingTwinChoices(*result.space_group, sg_search); + std::optional model_sibling; + if (!sibling_choices.empty() && config_.model && !model_predicted && model_group) { + const auto siblings = SiblingTwinChoices( + *result.space_group, sg_search, + ModelRotationsOnDataAxes(*model_group, model_cell, gemmi::UnitCell(*result.consensus_cell))); + if (!siblings.empty()) + model_sibling = siblings.front(); + } const bool model_tests_twin = config_.model - && (twin_metric_subgroup.has_value() + && (twin_metric_subgroup.has_value() || model_sibling.has_value() || (!model_predicted && model_rotations > 0 && model_rotations < result.space_group->operations().sym_ops.size())); if (!geometry_prepass && result.space_group && sg_search.point_group_order > 1 @@ -3592,18 +3655,29 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { model_twin_choice = ChooseTwinSubgroup( *result.space_group, sg_search, std::nullopt, ModelRotationsOnDataAxes(*model_group, model_cell, gemmi::UnitCell(*result.consensus_cell))); + if (!model_twin_choice && model_sibling) { + model_twin_choice = model_sibling; + sibling_for_model = model_sibling->subgroup.xhm(); + } if (model_twin_choice && twin_choice && SortedRotations(model_twin_choice->subgroup.operations()) == SortedRotations(twin_choice->subgroup.operations())) model_twin_choice.reset(); + if (model_sibling && model_twin_choice && model_twin_choice->sibling) + logger.Info("Twin hypothesis: the model's group {} is not a subgroup of {} but another index-2 " + "subgroup of the refused {}: perfectly twinned by {}, it merges as {} does - the " + "merge in {} is made for the model to test", model_group->xhm(), + result.space_group->xhm(), sg_search.refused_point_group_hm, + model_sibling->law.as_hkl().triplet(), result.space_group->xhm(), + model_sibling->subgroup.xhm()); + else if (model_twin_choice) + logger.Info("Twin hypothesis: the model's group {} asks for the merge in {}, with twin law {}, " + "made for the model to test", model_group->xhm(), model_twin_choice->subgroup.xhm(), + model_twin_choice->law.as_hkl().triplet()); if (model_twin_choice && twin_for_mr_only) { twin_choice = std::move(model_twin_choice); model_twin_choice.reset(); } - if (model_twin_choice) - logger.Info("Twin hypothesis: the model's group {} asks for the merge in {}, with twin law {}, " - "made for the model to test", model_group->xhm(), model_twin_choice->subgroup.xhm(), - model_twin_choice->law.as_hkl().triplet()); } twin_choice_reason = twin_overruled_subgroup ? fmt::format("a twin gate refused it ({}) and was overruled", twin_overruled_reason) @@ -3852,6 +3926,17 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { twin_merge = twin_law_merge(*twin_choice); if (model_twin_choice) model_twin_merge = twin_law_merge(*model_twin_choice); + // A model without a group of its own to choose a sibling hypothesis by (a prediction, or coordinates without a + // crystal cell) is searched for by molecular replacement in every sibling's merge, beside the answer's. Each + // stays in the scaler's frame here; std::list, as the files made from them point into them. + std::list> sibling_merges; + const bool model_run_search_siblings = config_.model && write_output && !geometry_prepass && !superseded + && (model_predicted || (config_.force_mr && !model_group)); + if (model_run_search_siblings) + for (const auto &c : sibling_choices) + if (!(model_twin_choice && model_twin_choice->subgroup.xhm() == c.subgroup.xhm())) + if (auto m = twin_law_merge(c)) + sibling_merges.emplace_back(c, std::move(*m)); std::optional ledger_p1; Logger ledger_log = Logger::Buffered(); std::future ledger_ahead; @@ -4066,7 +4151,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // on a copy given the same amplitudes. Its log is held, and printed where it is taken. Logger model_search_log = Logger::Buffered(); std::future model_search_run; - if (model_run && (model_predicted || config_.force_mr) && !twin_merge && !model_twin_merge + if (model_run && (model_predicted || config_.force_mr) && !twin_merge && !model_twin_merge && sibling_merges.empty() && config_.reference_data.empty()) { model_merged = sm.merged; model_search_run = std::async(std::launch::async, @@ -4571,10 +4656,12 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { gemmi::Op law; TwinLawFraction fraction; const RotationScaleMerge::Result *merge = nullptr; // the merge it was made from, for its statistics + bool sibling = false; // TwinSubgroupChoice::sibling }; // The model's own subgroup's merge (model_twin_choice), made ready the same way, and reported in the // model's section only: written, as _mr, only where the model chooses it. std::optional twin_file, model_twin_file; + std::list sibling_files; // the sibling hypotheses' merges molecular replacement searches (sibling_merges) std::optional model_twin_report; const auto make_twin_file = [&](const TwinSubgroupChoice &choice, const RotationScaleMerge::Result &twin_merge, const std::string &file_suffix, @@ -4607,7 +4694,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e)) c.subgroup = *other; } - if (!rotations_in(c.subgroup, adopted)) { + if (!c.sibling && !rotations_in(c.subgroup, adopted)) { logger.Warning("Twin hypothesis: {} is not a subgroup of {} as written; its merge is not written", c.subgroup.xhm(), adopted.xhm()); return std::nullopt; @@ -4661,7 +4748,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { logger.Info("Twin hypothesis: {} in {}, twin law {} at a fraction of at least {:.2f} ({} pairs) -> {}.mtz/.cif{}", adopted.xhm(), c.subgroup.xhm(), t.twin_law, fraction.fraction, fraction.pairs, t.file_suffix, file_suffix.empty() ? "" : " where the model chooses it"); - return TwinFile{std::move(merged), cell, c.subgroup, c.law, fraction, &twin_merge}; + return TwinFile{std::move(merged), cell, c.subgroup, c.law, fraction, &twin_merge, c.sibling}; }; const auto prepare_twin_file = [&] { if (!twin_file && twin_merge && twin_choice) { @@ -4696,6 +4783,18 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { std::optional placed_model_out; // The merge that placement was made against, which refinement refines it against. const std::vector *placed_merged = nullptr; + // The other of the two hypotheses the model was scored on - the answer's merge and the twin hypothesis' - + // where it fits both within HYPOTHESIS_R_FREE_MARGIN: refined beside the chosen one, which it replaces + // where it refines better by REFINEMENT_PATH_NOISE in R-free. `prior` holds the twin hypothesis' files + // where the rival is that hypothesis (its merged reflections are what it is refined against); the + // answer's merge otherwise. + struct RivalPlacement { + PlacedModelOutput placed; + const std::vector *merged = nullptr; + std::optional prior; + double r_free_validation = NAN; + }; + std::optional rival; // The twin hypothesis the model is asked about: its own subgroup's merge where it has one, the data's // otherwise - and whether that merge was made only for the model, so is written only where it chooses it. std::optional &model_twin = model_twin_merge ? model_twin_file : twin_file; @@ -4882,27 +4981,47 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { placed_model_out = PlacedModelOutput{twin_validation.placed_model, mt.cell, mt.space_group}; placed_merged = &mt.merged; } - // The model placed there and fitting better there chooses the subgroup: its files, _mr. They - // carry the model's own group where it has the subgroup's rotations - the screws and the hand are the - // model's, as the coordinates beside them are, unless the anomalous map at its atoms is inverted - and - // none of the reflections that group extinguishes. - if (maps_from_twin && placed_in_twin) { - twin_chosen_by_model = true; + // The model placed there chooses the subgroup where it also fits better there: its files, _mr. + // They carry the model's own group where it has the subgroup's rotations - the screws and the hand are + // the model's, as the coordinates beside them are, unless the anomalous map at its atoms is inverted - + // and none of the reflections that group extinguishes. The model is refined in that group, against + // those reflections: a sibling's merge (SiblingTwinChoices) is labelled with the symmorphic group of + // its rotations only. + std::optional twin_prior; + if (placed_in_twin && twin_validation.placed_model) { gemmi::SpaceGroup label = mt.space_group; - if (const gemmi::SpaceGroup *own = twin_validation.placed_model - ? twin_validation.placed_model->find_spacegroup() : nullptr; + if (const gemmi::SpaceGroup *own = twin_validation.placed_model->find_spacegroup(); own != nullptr && !twin_validation.anomalous_hands_disagree && SortedRotations(own->operations()) == SortedRotations(mt.space_group.operations())) label = *own; - ModelPriorFile f{mt.merged, mt.cell, label, mt.law, mt.fraction, twin_validation.placed_model, mt.merge}; + twin_prior = ModelPriorFile{mt.merged, mt.cell, label, mt.law, mt.fraction, + twin_validation.placed_model, mt.merge}; if (label.xhm() != mt.space_group.xhm()) { const gemmi::GroupOps lops = label.operations(); - std::erase_if(f.merged, [&](const MergedReflection &r) { + std::erase_if(twin_prior->merged, [&](const MergedReflection &r) { return lops.is_systematically_absent({{r.h, r.k, r.l}}); }); } - model_prior_file = std::move(f); + } + // The two hypotheses the model fits within HYPOTHESIS_R_FREE_MARGIN of each other are both carried into + // refinement, and decided there (see the refinement below): the one not chosen here is the rival. + const bool both_fit = validation.model_fits && validation.placed_model && twin_prior; + const bool close = both_fit && std::abs(twin_r_free(twin_validation) - answer_r_free) < HYPOTHESIS_R_FREE_MARGIN; + if (maps_from_twin && twin_prior) { + twin_chosen_by_model = true; + model_prior_file = std::move(twin_prior); result.model_prior_file_suffix = "_mr"; + if (placed_model_out) { + placed_model_out->space_group = model_prior_file->space_group; + placed_merged = &model_prior_file->merged; + } + if (close) + rival = RivalPlacement{PlacedModelOutput{validation.placed_model, *result.consensus_cell, + experiment_.GetSpaceGroupOrP1()}, + &sm.merged, std::nullopt, answer_r_free}; + } else if (close && !maps_from_twin) { + PlacedModelOutput p{twin_prior->placed, twin_prior->cell, twin_prior->space_group}; + rival = RivalPlacement{std::move(p), nullptr, std::move(twin_prior), twin_r_free(twin_validation)}; } mt_report.maps_from_this_merge = maps_from_twin; mt_report.model_validation = std::move(twin_validation); @@ -4932,6 +5051,8 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { rsm.reset(); const size_t merge_on_card = rsm && rsm->OnGPU() ? rsm->DeviceBytes() : 0; MolecularReplacementResult mr; + // The merge each candidate's `merge` index names, past the answer's (index 0, null here). + std::vector mr_twin_files{nullptr}; if (model_search_run.valid()) { TimingMark wait_mark(logger, "molecular replacement (beside the analyses) wait"); mr = model_search_run.get(); @@ -4939,9 +5060,21 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { model_search_log.ReplayInto(logger); } else { std::vector merges{{"MERGE", &sm.merged, *result.consensus_cell, experiment_.GetSpaceGroupOrP1()}}; - if (model_twin) + if (model_twin) { merges.push_back({"TWIN_SUBGROUP_MERGE", &model_twin->merged, model_twin->cell, - model_twin->space_group, model_twin->law}); + model_twin->space_group, model_twin->law, model_twin->sibling}); + mr_twin_files.push_back(&*model_twin); + } + // Every sibling hypothesis, for a model without a group to choose one by (sibling_merges). + for (auto &[c, m] : sibling_merges) { + ProcessResult::TwinSubgroupMerge t; + if (auto f = make_twin_file(c, m, "_mr", t)) + sibling_files.push_back(std::move(*f)); + } + for (const TwinFile &f : sibling_files) { + merges.push_back({"SIBLING_MERGE", &f.merged, f.cell, f.space_group, f.law, true}); + mr_twin_files.push_back(&f); + } mr = search_model(merges, logger); } if (mr.chosen >= 0) { @@ -4969,10 +5102,18 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { return gops.is_systematically_absent({{r.h, r.k, r.l}}); }); model_prior_file = std::move(f); - } else if (mr.decision == "TWIN_SUBGROUP" && model_twin) { - model_prior_file = ModelPriorFile{model_twin->merged, model_twin->cell, c.space_group, model_twin->law, - model_twin->fraction, c.placed, model_twin->merge}; - twin_chosen_by_model = true; + } else if (mr.decision == "TWIN_SUBGROUP" && c.merge > 0 + && static_cast(c.merge) < mr_twin_files.size()) { + const TwinFile &tf = *mr_twin_files[c.merge]; + model_prior_file = ModelPriorFile{tf.merged, tf.cell, c.space_group, tf.law, tf.fraction, c.placed, + tf.merge}; + if (c.space_group.xhm() != tf.space_group.xhm()) { + const gemmi::GroupOps gops = c.space_group.operations(); + std::erase_if(model_prior_file->merged, [&](const MergedReflection &r) { + return gops.is_systematically_absent({{r.h, r.k, r.l}}); + }); + } + twin_chosen_by_model = model_twin && &tf == &*model_twin; } result.model_prior_file_suffix = model_prior_file ? "_mr" : ""; // The maps, the R-factors and CC(model, data): the placement validated against the merge it was @@ -5001,9 +5142,12 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { v.r_work_sigma = rb.r_work_sigma; v.null_replicates = rb.null_replicates; } - if (own) + // Refined where it was chosen: in the group of the model's choice, against that choice's merge. + if (own) { model_prior_file->placed = v.placed_model; - else { + placed_model_out = PlacedModelOutput{v.placed_model, cell, data_sg}; + placed_merged = &model_prior_file->merged; + } else { placed_model_out = PlacedModelOutput{v.placed_model, cell, data_sg}; placed_merged = &sm.merged; } @@ -5043,6 +5187,41 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { }); } + // The symmetry hypotheses a refused promotion leaves (SiblingTwinChoices), on the axes the answer is written + // on: report only, whether or not a model tested one. + if (!sibling_choices.empty() && result.space_group_search + && result.space_group_search->refused_point_group_representative) { + const auto moved = [&](TwinSubgroupChoice c) { + if (!(merge_to_written == gemmi::Op::identity())) + MoveTwinSubgroupChoice(c, merge_to_written.inverse()); + return c; + }; + const auto &search = *result.space_group_search; + ProcessResult::SymmetryHypothesis adopted{experiment_.GetSpaceGroupOrP1().xhm(), {}, "ADOPTED", false}; + adopted.tested_by_model = result.model_validation && result.model_validation->ok; + result.symmetry_hypotheses.push_back(adopted); + // The index-2 subgroup of the refused group that is not a sibling has the adopted group's rotations: it, + // twinned by a law of the refused group outside it. + for (const auto &t : IndexTwoTwinSubgroups(*search.refused_point_group_representative)) + if (std::none_of(sibling_choices.begin(), sibling_choices.end(), [&](const TwinSubgroupChoice &c) { + return SortedRotations(c.subgroup.operations()) == SortedRotations(t.subgroup.operations()); + })) { + const TwinSubgroupChoice c = moved({t.subgroup, t.laws.front(), {}, false}); + result.symmetry_hypotheses.push_back( + {adopted.space_group, c.law.as_hkl().triplet('h'), "ADOPTED_TWINNED", adopted.tested_by_model}); + } + for (const auto &sc : sibling_choices) { + const TwinSubgroupChoice c = moved(sc); + const bool tested = (!sibling_files.empty() && result.model_prior) + || sibling_for_model == sc.subgroup.xhm() + && ((model_twin_report && model_twin_report->model_validation) + || (result.twin_subgroup_merge && result.twin_subgroup_merge->model_validation) + || result.model_prior); + result.symmetry_hypotheses.push_back({c.subgroup.xhm(), c.law.as_hkl().triplet('h'), "SIBLING_TWINNED", + tested}); + } + } + if (model_run) { result.model_placement = placed_by_rigid_body ? "RIGID_BODY" : result.model_prior && result.model_prior->chosen >= 0 ? "MR" : "NOT_PLACED"; @@ -5054,15 +5233,6 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { why = "the model could not be searched for: " + result.model_prior->failure_reason; if (!(result.model_validation && result.model_validation->ok) && !why.empty()) result.warnings.push_back({PathologyCode::MODEL_NOT_VALIDATED, "Model validation did not run: " + why}); - // The twin hypothesis' merge made only for the model to test: written, and reported, only where it chose it. - if (model_twin_only && !twin_chosen_by_model) { - model_twin.reset(); - model_twin_entry.reset(); - } - result.model_twin_subgroup_merge = model_twin_report; - if (model_prior_file) - result.model_prior_file_space_group = model_prior_file->space_group.xhm(); - result.outputs->placed_model = placed_model_out; // The outcome, as notes of the step that decided it if (result.model_prior) { const MolecularReplacementResult &mp = *result.model_prior; @@ -5100,7 +5270,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { settings.gpu_resident_bytes = rsm && rsm->OnGPU() ? rsm->DeviceBytes() : 0; // 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; + const PlacedModelOutput p = *placed_model_out; // processing-inert: a refinement that fails is reported, and the run goes on to its other outputs RefinementResult refined; try { @@ -5110,8 +5280,66 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { refined = RefinementResult(); refined.failure_reason = std::string("refinement failed: ") + e.what(); } + // The rival hypothesis (see RivalPlacement), refined the same way, and the two decided on the refined + // R-free: each refinement has made its own twin decision against its own null, so what is compared is + // the best each hypothesis can do. The validation's choice stands unless the rival beats it by + // REFINEMENT_PATH_NOISE. + if (rival && refined.ok && !refined.cancelled) { + phase("Refining the model (rival symmetry hypothesis)"); + const RivalPlacement &rv = *rival; + const std::vector &rv_merged = rv.prior ? rv.prior->merged : *rv.merged; + RefinementResult other; + try { + other = RefineModel(*rv.placed.model, rv_merged, gemmi::UnitCell(rv.placed.cell), + rv.placed.space_group, experiment_.GetWavelength_A(), logger, settings); + } catch (const std::exception &e) { + other = RefinementResult(); + other.failure_reason = std::string("refinement failed: ") + e.what(); + } + ProcessResult::RivalRefinement rr; + rr.space_group = rv.placed.space_group.xhm(); + rr.merge = rv.prior ? "TWIN_SUBGROUP_MERGE" : "MERGE"; + rr.chosen_space_group = p.space_group.xhm(); + rr.chosen_r_free = refined.r_free; + if (other.ok) { + rr.r_work = other.r_work; + rr.r_free = other.r_free; + rr.twin = other.twin_decision; + rr.twin_law = other.twin_law; + } + rr.replaced = other.ok && !other.cancelled && other.r_free < refined.r_free - REFINEMENT_PATH_NOISE; + logger.Info("Symmetry hypotheses: refined in {} R-free {:.4f}{}, in {} (the rival) {} - {}", + p.space_group.xhm(), refined.r_free, + refined.twin_decision == "TWINNED" ? " twinned by " + refined.twin_law : std::string(), + rr.space_group, + other.ok ? fmt::format("R-free {:.4f}{}", other.r_free, + other.twin_decision == "TWINNED" ? " twinned by " + other.twin_law + : std::string()) + : other.failure_reason, + rr.replaced ? "the rival is taken" : "the validation's choice stands"); + if (rr.replaced) { + // The model's choice moves with it: the twin hypothesis' files where the rival is it, none + // (the answer's merge) where it is the answer. + if (rival->prior) { + model_prior_file = std::move(rival->prior); + result.model_prior_file_suffix = "_mr"; + twin_chosen_by_model = true; + placed_model_out = PlacedModelOutput{model_prior_file->placed, model_prior_file->cell, + model_prior_file->space_group}; + } else { + model_prior_file.reset(); + result.model_prior_file_suffix.clear(); + twin_chosen_by_model = false; + placed_model_out = rival->placed; + } + refined = std::move(other); + } + result.rival_refinement = std::move(rr); + result.model_choice_decided_by = "REFINEMENT"; + } + const PlacedModelOutput &p_final = *placed_model_out; if (refined.ok) { - result.outputs->refined_model = PlacedModelOutput{refined.model, p.cell, p.space_group}; + result.outputs->refined_model = PlacedModelOutput{refined.model, p_final.cell, p_final.space_group}; result.outputs->refined_maps = refined.maps; result.outputs->blobs = refined.blobs; result.outputs->solvent = refined.solvent; @@ -5123,6 +5351,20 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { logger.Warning("Refinement did not run: {}", refined.failure_reason); result.refinement = std::move(refined); } + // The twin hypothesis' merge made only for the model to test: written, and reported, only where it chose it. + if (model_twin_only && !twin_chosen_by_model) { + model_twin.reset(); + model_twin_entry.reset(); + } + result.model_twin_subgroup_merge = model_twin_report; + if (model_prior_file) { + result.model_prior_file_space_group = model_prior_file->space_group.xhm(); + if (model_prior_file->law) + result.model_prior_file_twin_law = model_prior_file->law->as_hkl().triplet('h'); + } + if (result.model_choice_decided_by.empty() && result.model_validation && result.model_validation->ok) + result.model_choice_decided_by = "VALIDATION"; + result.outputs->placed_model = placed_model_out; } if (cell_esd.valid()) { diff --git a/tests/TwinningAnalysisTest.cpp b/tests/TwinningAnalysisTest.cpp index fb0a48862..f47fd915f 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -476,6 +476,40 @@ TEST_CASE("Index-2 subgroups with a two-fold twin law", "[twinning]") { CHECK(names("C 2 2 21").size() == 3); } +// The perfect-twin hypotheses a refused promotion leaves beside the adopted group: the other index-2 +// subgroups of the refused group, each with the two-folds of its coset as laws. A P3x21 adopted under a +// refused 622 leaves P6 (law k,h,-l among others) and P312; a group that is not an index-2 subgroup of +// the refused one leaves none. +TEST_CASE("Sibling twin hypotheses of a refused promotion", "[twinning]") { + const auto siblings = [](const char *adopted, const char *refused) { + return SiblingTwinSubgroups(gemmi::get_spacegroup_by_name(adopted), gemmi::get_spacegroup_by_name(refused)); + }; + const auto names = [](const std::vector &v) { + std::vector out; + for (const auto &t : v) { + CHECK_FALSE(t.laws.empty()); + out.push_back(t.subgroup.xhm()); + } + std::sort(out.begin(), out.end()); + return out; + }; + const auto hexagonal = siblings("P 31 2 1", "P 6 2 2"); + CHECK(names(hexagonal) == std::vector{"P 3 1 2", "P 6"}); + for (const auto &t : hexagonal) { + std::vector laws; + for (const auto &l : t.laws) + laws.push_back(l.as_hkl().triplet('h')); + CAPTURE(t.subgroup.xhm(), laws); + if (t.subgroup.xhm() == "P 6") + CHECK(std::find(laws.begin(), laws.end(), "k,h,-l") != laws.end()); + else + CHECK(std::find(laws.begin(), laws.end(), "-h,-k,l") != laws.end()); + } + CHECK(names(siblings("P 1 21 1", "P 21 21 21")) == std::vector{"P 1 1 21", "P 21 1 1"}); + CHECK(siblings("P 4", "P 6 2 2").empty()); + CHECK(siblings("P 6 2 2", "P 6 2 2").empty()); +} + // Yeates' H on the pairs a law relates, in the group the law is not part of, reads the twin fraction. TEST_CASE("Twin fraction of a law from its pairs' disagreement", "[twinning]") { const gemmi::SpaceGroup &p4 = gemmi::get_spacegroup_by_name("P 4"); From 821636f2939d95d49e74a593f64b66e312b71239 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 13:35:41 +0100 Subject: [PATCH 2/3] RESULT-CHANGING: the placed model is refined as the merge's group holds it; rivals only where both hold it whole A placement was refined in the merge's label group. Where the label leaves a screw open and the model asserts one with the same rotations (a P6 label with a P63 model, P321 with P3121), the copies were expanded by the wrong operators; where the model's group has more rotations than a twin hypothesis' subgroup, only part of the crystal's content was refined. RefinableIn: the model's own group where it has the merge's rotations (unless the anomalous map is inverted), the model expanded into the merge's group (ExpandIntoSubgroup) where its group is larger. The best-two rival refinement is made only where both merges hold the model whole that way: on the subset battery a P63 model refined in a P6 label lost to its own P3 expansion-less twin rival before this. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- docs/RUGNUX_REPORT.md | 8 +++- rugnux/RugnuxScaleMerge.cpp | 73 ++++++++++++++++++++++++++++++------- 2 files changed, 65 insertions(+), 16 deletions(-) diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 4820546b9..4d7e759d1 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -1009,8 +1009,12 @@ the model is scored on, exactly as a subgroup's above, and the model is refined With a model that has no group of its own (a prediction, or `--mr` with coordinates without a cell), every sibling's merge is searched by molecular replacement beside the answer's. -Where the model fits both the answer's merge and the twin hypothesis' within 0.03 in R-free before -refinement, **both are refined**, each with its own twin decision (below), and the refined R-free decides +The placed model is refined as the merge's group holds it: in the model's own group where that has the +merge's rotations (so the screws the merge's label leaves open are the model's), and expanded into the +merge's group where the model's group has more rotations (each symmetry image the merge's group lacks +becomes an independent copy). Where the model fits both the answer's merge and the twin hypothesis' within +0.03 in R-free before refinement, and each merge's group holds the model whole that way, **both are +refined**, each with its own twin decision (below), and the refined R-free decides which the model chooses: the rival replaces the validation's choice only where it refines better by 0.01, what two refinements of one model differ by from their paths alone. `MODEL_CHOICE_DECIDED_BY=` is `VALIDATION` or `REFINEMENT`; the rival's `RIVAL_REFINE_SPACE_GROUP`, `RIVAL_REFINE_MERGE`, diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 78eec1738..9746b96f0 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -422,6 +422,42 @@ namespace { return out; } + // A placement as it is refined against a merge in `merge_sg`, cell `cell`: in the model's own group where that + // has the merge's rotations - the model's screws, which the merge's label may leave open, unless the anomalous + // map at its atoms is inverted; expanded into the merge's group where the model's group has more rotations, + // every symmetry image the merge's group lacks an independent copy (ExpandIntoSubgroup), as a crystal of the + // merge's group; as it is otherwise, which the merge's group holds only in part (`holds` false). + struct RefinablePlacement { + std::shared_ptr model; + gemmi::SpaceGroup space_group; + bool holds = false; + }; + RefinablePlacement RefinableIn(const std::shared_ptr &placed, const UnitCell &cell, + const gemmi::SpaceGroup &merge_sg, bool hands_disagree) { + RefinablePlacement out{placed, merge_sg, false}; + const gemmi::SpaceGroup *own = placed ? placed->find_spacegroup() : nullptr; + if (own == nullptr) + return out; + const auto m = SortedRotations(merge_sg.operations()); + const auto o = SortedRotations(own->operations()); + if (o == m) { + if (!hands_disagree) + out.space_group = *own; + out.holds = true; + } else if (o.size() > m.size() && std::includes(o.begin(), o.end(), m.begin(), m.end())) { + const gemmi::UnitCell h(cell); + if (const auto setting = SubgroupBasis(h, *own, h, merge_sg)) { + int n = 0; + out.model = std::make_shared( + ExpandIntoSubgroup(*placed, 1, *own, *setting, h, merge_sg, &n)); + out.holds = n > 1; + if (!out.holds) + out.model = placed; + } + } + return out; + } + // A twin-law choice on other axes: the group by SpaceGroupInBasis, the law conjugated by the // Miller-index operator the change of basis moves the reflections by. False where no setting in // GEMMI's table names the moved group. @@ -4988,16 +5024,14 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // those reflections: a sibling's merge (SiblingTwinChoices) is labelled with the symmorphic group of // its rotations only. std::optional twin_prior; + bool twin_holds = false; if (placed_in_twin && twin_validation.placed_model) { - gemmi::SpaceGroup label = mt.space_group; - if (const gemmi::SpaceGroup *own = twin_validation.placed_model->find_spacegroup(); - own != nullptr && !twin_validation.anomalous_hands_disagree - && SortedRotations(own->operations()) == SortedRotations(mt.space_group.operations())) - label = *own; - twin_prior = ModelPriorFile{mt.merged, mt.cell, label, mt.law, mt.fraction, - twin_validation.placed_model, mt.merge}; - if (label.xhm() != mt.space_group.xhm()) { - const gemmi::GroupOps lops = label.operations(); + const RefinablePlacement r = RefinableIn(twin_validation.placed_model, mt.cell, mt.space_group, + twin_validation.anomalous_hands_disagree); + twin_holds = r.holds; + twin_prior = ModelPriorFile{mt.merged, mt.cell, r.space_group, mt.law, mt.fraction, r.model, mt.merge}; + if (r.space_group.xhm() != mt.space_group.xhm()) { + const gemmi::GroupOps lops = r.space_group.operations(); std::erase_if(twin_prior->merged, [&](const MergedReflection &r) { return lops.is_systematically_absent({{r.h, r.k, r.l}}); }); @@ -5005,19 +5039,26 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { } // The two hypotheses the model fits within HYPOTHESIS_R_FREE_MARGIN of each other are both carried into // refinement, and decided there (see the refinement below): the one not chosen here is the rival. - const bool both_fit = validation.model_fits && validation.placed_model && twin_prior; + // Only where each merge's group holds the model whole (RefinableIn): a refinement of part of the + // crystal's content is no measure of its hypothesis. + const RefinablePlacement answer_refinable = + RefinableIn(validation.placed_model, *result.consensus_cell, experiment_.GetSpaceGroupOrP1(), + validation.anomalous_hands_disagree); + const bool both_fit = validation.model_fits && validation.placed_model && twin_prior + && twin_holds && answer_refinable.holds; const bool close = both_fit && std::abs(twin_r_free(twin_validation) - answer_r_free) < HYPOTHESIS_R_FREE_MARGIN; if (maps_from_twin && twin_prior) { twin_chosen_by_model = true; model_prior_file = std::move(twin_prior); result.model_prior_file_suffix = "_mr"; if (placed_model_out) { + placed_model_out->model = model_prior_file->placed; placed_model_out->space_group = model_prior_file->space_group; placed_merged = &model_prior_file->merged; } if (close) - rival = RivalPlacement{PlacedModelOutput{validation.placed_model, *result.consensus_cell, - experiment_.GetSpaceGroupOrP1()}, + rival = RivalPlacement{PlacedModelOutput{answer_refinable.model, *result.consensus_cell, + answer_refinable.space_group}, &sm.merged, std::nullopt, answer_r_free}; } else if (close && !maps_from_twin) { PlacedModelOutput p{twin_prior->placed, twin_prior->cell, twin_prior->space_group}; @@ -5026,9 +5067,13 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { mt_report.maps_from_this_merge = maps_from_twin; mt_report.model_validation = std::move(twin_validation); } + // Refined as the answer's merge holds it (RefinableIn): in the model's own group where it has the merge's + // rotations, so that screws the merge's label leaves open are the model's. if (validation.placed_model && !maps_from_twin) { - placed_model_out = PlacedModelOutput{validation.placed_model, *result.consensus_cell, - experiment_.GetSpaceGroupOrP1()}; + const RefinablePlacement r = RefinableIn(validation.placed_model, *result.consensus_cell, + experiment_.GetSpaceGroupOrP1(), + validation.anomalous_hands_disagree); + placed_model_out = PlacedModelOutput{r.model, *result.consensus_cell, r.space_group}; placed_merged = &sm.merged; } } From 18263bccac159a86e38aee8dd9ae196966e09e17 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 14:10:15 +0100 Subject: [PATCH 3/3] SiblingTwinSubgroups: a centred crystal's siblings carry its centring The refused group's representative can be given on the primitive lattice, which named primitive siblings (P222 beside an adopted I41). The refused group is taken with the adopted group's centring before its index-2 subgroups are listed. Report-only on every set the subset battery holds; a sibling-group model of a centred crystal is now matched. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- image_analysis/scale_merge/TwinningAnalysis.cpp | 7 ++++++- tests/TwinningAnalysisTest.cpp | 3 +++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 79af2e7fb..fe7a27109 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -72,8 +72,13 @@ std::vector SiblingTwinSubgroups(const gemmi::SpaceGroup& adopted, const auto g = rotations(refused); if (2 * a.size() != g.size() || !std::includes(g.begin(), g.end(), a.begin(), a.end())) return {}; + // The refused group with the adopted one's centring: a representative given on the primitive lattice + // would name primitive siblings of a centred crystal. + gemmi::GroupOps centred = refused.operations(); + centred.cen_ops = adopted.operations().cen_ops; + const gemmi::SpaceGroup* parent = gemmi::find_spacegroup_by_ops(centred); std::vector out; - for (auto& t : IndexTwoTwinSubgroups(refused)) + for (auto& t : IndexTwoTwinSubgroups(parent != nullptr ? *parent : refused)) if (rotations(t.subgroup) != a) out.push_back(std::move(t)); return out; diff --git a/tests/TwinningAnalysisTest.cpp b/tests/TwinningAnalysisTest.cpp index f47fd915f..49141357b 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -506,6 +506,9 @@ TEST_CASE("Sibling twin hypotheses of a refused promotion", "[twinning]") { CHECK(std::find(laws.begin(), laws.end(), "-h,-k,l") != laws.end()); } CHECK(names(siblings("P 1 21 1", "P 21 21 21")) == std::vector{"P 1 1 21", "P 21 1 1"}); + // A centred crystal's siblings are centred, whatever centring the refused group was given with. + const auto centred = names(siblings("I 41", "P 4 2 2")); + CHECK(std::find(centred.begin(), centred.end(), "I 2 2 2") != centred.end()); CHECK(siblings("P 4", "P 6 2 2").empty()); CHECK(siblings("P 6 2 2", "P 6 2 2").empty()); }