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 9dba7759b..e3b9ceb6f 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -997,6 +997,33 @@ 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. + +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`, +`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..fe7a27109 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -60,6 +60,30 @@ 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 {}; + // 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(parent != nullptr ? *parent : 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 bed0cc371..1a85a59ba 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)); } @@ -2744,6 +2763,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 42ee54ccd..f240ce0f3 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -759,6 +759,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 685d48615..0d7d38991 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,72 @@ 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 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. @@ -1996,6 +2078,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. @@ -3581,9 +3667,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 @@ -3601,18 +3700,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) @@ -3861,6 +3971,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; @@ -4075,7 +4196,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, @@ -4580,10 +4701,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, @@ -4616,7 +4739,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; @@ -4670,7 +4793,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) { @@ -4705,6 +4828,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; @@ -4936,34 +5071,63 @@ 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; - gemmi::SpaceGroup label = mt.space_group; - if (const gemmi::SpaceGroup *own = twin_validation.placed_model - ? twin_validation.placed_model->find_spacegroup() : nullptr; - 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}; - if (label.xhm() != mt.space_group.xhm()) { - const gemmi::GroupOps lops = label.operations(); - std::erase_if(f.merged, [&](const MergedReflection &r) { + // 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; + bool twin_holds = false; + if (placed_in_twin && twin_validation.placed_model) { + 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}}); }); } - 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. + // 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{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}; + 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); } + // 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; } } @@ -4986,6 +5150,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(); @@ -4993,9 +5159,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) { @@ -5024,10 +5202,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 @@ -5056,9 +5242,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; } @@ -5108,6 +5297,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"; @@ -5119,15 +5343,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; @@ -5165,22 +5380,20 @@ 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; - // 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; + const PlacedModelOutput p = *placed_model_out; // a copy: a rival taken below replaces *placed_model_out + // In the group RefinableIn chose: the model's own where it has the rotations of the group written - + // the screws are then the model's, as its fit already took them - so 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), does not pack its copies + // wrongly. The reflections that group extinguishes are left out. + const 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())) { + if (const gemmi::SpaceGroup written = experiment_.GetSpaceGroupOrP1(); + refine_sg.xhm() != written.xhm() + && SortedRotations(refine_sg.operations()) == SortedRotations(written.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(); + refine_sg.xhm(), written.xhm()); + const gemmi::GroupOps ops = refine_sg.operations(); for (const MergedReflection &r : *placed_merged) if (!ops.is_systematically_absent({{r.h, r.k, r.l}})) refine_merged.push_back(r); @@ -5195,8 +5408,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, refine_sg}; + 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; @@ -5208,6 +5479,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..49141357b 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -476,6 +476,43 @@ 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"}); + // 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()); +} + // 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");