From 097668e674e3b21e5f8e324e87e546dc5c67b8c6 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 17:12:55 +0200 Subject: [PATCH 1/4] rugnux: keep the second pass's supercell when its own axis-harmonic arbiter settled it The two-pass supercell guards (in the first-pass indexing and after the second pass) revert to pass 1's lattice whenever the second pass's primitive cell is markedly larger. On a crystal with a doubled axis pass 1's schemes did not tie on frames, so pass 1 never asked the axis-harmonic arbiter and kept the halved cell; pass 2's schemes tied, its arbiter established the larger cell on the spots (17048 against 16099 validation spots), and the guard then threw it away. The scheme comparison now records when the lattice it keeps is the larger cell of an axis-harmonic pair the arbiter settled, and both guards leave such a lattice in place. Where the arbiter did not settle the larger cell - every other set the guard fires on - nothing changes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- rugnux/Rugnux.h | 3 +++ rugnux/RugnuxFirstPass.cpp | 3 +++ rugnux/RugnuxFirstPass.h | 4 +++- rugnux/RugnuxFirstPassAccept.cpp | 11 ++++++++++- rugnux/RugnuxPasses.cpp | 2 +- 5 files changed, 20 insertions(+), 3 deletions(-) diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 1147a40f1..4b5c3dd46 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -684,6 +684,9 @@ class Rugnux { // force_rotation_result_), so it integrates on the correct cell with pass-1's refined orientation/axis - a // lattice-only force loses that metadata and integrates markedly worse. std::optional prepass_result_; + // The last second pass found a supercell of pass 1's lattice and its own axis-harmonic arbiter + // established that larger cell on the spots, so neither supercell guard reverts to pass 1. + bool supercell_settled_ = false; // What the pre-scan measured about powder contamination - a crystalline phase other than the // crystal, diffracting as rings. Measured on every run and reported whatever it says; acted on // only where a pass that failed indexes better with the rings set aside. diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index dade24cc5..4d52b6613 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -974,6 +974,7 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin name, vol, bp.name, bp.vol, harmonic_ratio); if (!bp.result.has_value() || clearly_more || integer_subcell) { + bp.larger_settled = harmonic_settled && vol > bp.vol; if (integer_subcell && vol < bp.vol) logger.Info("Scheme '{}' cell (vol {:.0f}) is a {:.0f}x sub-cell of '{}' (vol {:.0f}) at " "equal frame count - adopting the smaller true cell", name, vol, @@ -986,6 +987,8 @@ auto Rugnux::FirstPassRun::ScoreSchemes(const SchemeIndexers &ris, IndexAndRefin bp.vol = vol; bp.name = name; bp.result = std::move(result); + } else if (harmonic_settled) { + bp.larger_settled = true; } } return bp; diff --git a/rugnux/RugnuxFirstPass.h b/rugnux/RugnuxFirstPass.h index a268c97eb..6909e172d 100644 --- a/rugnux/RugnuxFirstPass.h +++ b/rugnux/RugnuxFirstPass.h @@ -28,8 +28,10 @@ namespace rugnux_internal { // A candidate lattice of the first pass and what it scored on the validation frames. // indexer_error: what the indexer threw, if it threw. An empty result with no error is a // verdict on the crystal; an empty result with one is a verdict on the machine. + // larger_settled: the result is the larger cell of an axis-harmonic pair, and the arbiter in + // the scheme comparison established it on the spots. struct FirstPass { std::optional result; int score = -1; std::string name; double vol = 0.0; - std::optional indexer_error; }; + std::optional indexer_error; bool larger_settled = false; }; } struct Rugnux::FirstPassRun { diff --git a/rugnux/RugnuxFirstPassAccept.cpp b/rugnux/RugnuxFirstPassAccept.cpp index abca5c1bf..02bc68706 100644 --- a/rugnux/RugnuxFirstPassAccept.cpp +++ b/rugnux/RugnuxFirstPassAccept.cpp @@ -487,7 +487,16 @@ void Rugnux::FirstPassRun::Pass1SupercellGuard(const FirstPass &best) { const bool lost_centring = rx.prepass_result_->search_result.centering != 'P' && best.result->search_result.centering == 'P' && v1 > 1.0 && std::abs(v2 / v1 - 1.0) < 0.1; - if (supercell || lost_centring) { + // ... unless this pass's own axis-harmonic arbiter settled the larger cell on the spots. + // Pass 1 then never asked that question - its schemes did not tie on frames, so the + // arbiter did not run - and its smaller cell is the axis harmonic, not the reference + // (measured: a doubled axis the deposition confirms, thrown away here). + rx.supercell_settled_ = supercell && best.larger_settled; + if (rx.supercell_settled_) + logger.Info("Two-pass: that lattice ({:.0f} A^3) is {:.2f}x pass 1's ({:.0f} A^3), and " + "this pass's axis-harmonic arbiter established the larger cell on the " + "spots - keeping it", v2, v2 / v1, v1); + else if (supercell || lost_centring) { logger.Info("Two-pass: that lattice ({}-centred, volume {:.0f} A^3) disagrees with " "pass 1 ({}-centred, {:.0f} A^3) - integrating with pass-1's lattice instead", best.result->search_result.centering, v2, diff --git a/rugnux/RugnuxPasses.cpp b/rugnux/RugnuxPasses.cpp index ad420551e..ae893da22 100644 --- a/rugnux/RugnuxPasses.cpp +++ b/rugnux/RugnuxPasses.cpp @@ -1234,7 +1234,7 @@ std::optional Rugnux::RefinedPasses(RugnuxObserver *observer, Pro }; const double v1 = primitive_volume(*pass1.consensus_cell, pass1.consensus_centering); const double v2 = primitive_volume(*pass2.consensus_cell, pass2.consensus_centering); - if (v1 > 0.0 && v2 > 1.5 * v1) { + if (v1 > 0.0 && v2 > 1.5 * v1 && !supercell_settled_) { logger.Info("Two-pass: second-pass primitive cell volume {:.0f} A^3 is {:.2f}x pass-1's " "({:.0f} A^3) - a spurious supercell; re-running the second pass with " "pass-1's result forced", From 8f05793b670ef52ffcc6384bab35b7ba7fd48a02 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 17:15:00 +0200 Subject: [PATCH 2/4] Twinning report: decide whether a twin law is possible from the lattice Whether a twin law can exist was read off a hard-coded list of Laue classes (4/mmm, 6/mmm, m-3m and rhombohedral -3m "holohedral", everything else eligible). That is a property of the lattice, not of the class: a C-centred orthorhombic or a monoclinic lattice with no special metric admits no twin law for mmm or 2/m either, and a pseudo-merohedral one admits a law the class says nothing about. On five pseudo-merohedral sets the report said "Twinning INDICATED - refinement needs a twin law" where the lattice admits none. TwinLawPossible asks gemmi's find_twin_laws (Le Page two-folds within 3 deg, as RfreeFlags and ReindexAmbiguity already do), and AnalyzeTwinning takes the cell to decide it; in P1 or without a cell a law stays possible, so the pre-promotion analysis and the space-group search are unchanged. Twin-like statistics (<|L|> below 0.44) above the false-operator bar where no law is possible are now reported as NO_TWIN_LAW: not a twin in this group, so the symmetry may be too high or something else narrowed the distribution. Report-only: no reflection file changes. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- .../scale_merge/TwinningAnalysis.cpp | 62 ++++++++++--------- image_analysis/scale_merge/TwinningAnalysis.h | 27 ++++++-- rugnux/ResultReport.cpp | 27 ++++++-- rugnux/RugnuxScaleMerge.cpp | 5 +- rugnux/rugnux_cli.cpp | 6 +- tests/TwinningAnalysisTest.cpp | 53 +++++++++++++++- 6 files changed, 135 insertions(+), 45 deletions(-) diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index b8b07525c..01349f9a8 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -16,6 +16,8 @@ #include +#include "gemmi/twin.hpp" + #include "../../common/ParallelFor.h" namespace { @@ -28,12 +30,6 @@ namespace { return std::isfinite(r.I) && std::isfinite(r.d) && r.d > 0.0; } - // Merohedral twinning needs a twin law - a lattice symmetry operation that is not a symmetry of the - // crystal - which exists only when the Laue class is a proper subgroup of the lattice holohedry. - // The holohedral high-symmetry Laue classes (4/mmm, 6/mmm, m-3m, and -3m on a rhombohedral lattice) - // admit no such operation, so twinning is geometrically impossible and the intensity statistics - // cannot be indicating it. Low-symmetry classes stay eligible because pseudo-merohedral twinning - // through an accidental metric specialisation cannot be excluded from the symmetry alone. // <|L|> of a partial twin with fraction a: I = (1-a) X + a Y for independent Wilson X, Y, and the // expectation of |I1 - I2|/(I1 + I2) over two independent such intensities, integrated in closed // form (0.500 at a = 0, 0.440 at 0.1, 0.408 at 0.2, 0.375 at 0.5; checked by simulation, and @@ -60,26 +56,16 @@ namespace { } return 0.5 * (lo + hi); } +} - bool MerohedralTwinningPossible(const gemmi::SpaceGroup* sg) { - if (!sg) - return true; // P1 / unknown symmetry: cannot rule twinning out - switch (sg->laue_class()) { - case gemmi::Laue::L4mmm: // 4/mmm - tetragonal holohedry - case gemmi::Laue::L6mmm: // 6/mmm - hexagonal holohedry - case gemmi::Laue::Lm3m: // m-3m - cubic holohedry - return false; - case gemmi::Laue::L3m: // -3m is holohedral on a rhombohedral (R) lattice, but a - return sg->hm[0] != 'R'; // hexagonal-P 32/3m crystal can still twin towards 6/mmm - default: - return true; - } - } +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(); } TwinningAnalysisResult AnalyzeTwinning(const std::vector& merged, const gemmi::SpaceGroup* space_group, int resolution_shells, - const std::array* tncs_vector, char lattice_centring) { + const std::array* tncs_vector, char lattice_centring, + const gemmi::UnitCell* cell) { TwinningAnalysisResult result; if (merged.empty()) return result; @@ -326,7 +312,13 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg // ~490 A axis, whose neighbouring reflections overlap and share counts - averaging of another // kind, which narrows the distribution the same way. That is why this is a warning and never a // veto on the space group. - result.merohedral_twinning_possible = MerohedralTwinningPossible(space_group); + // Whether a twin law can exist is a property of the LATTICE: a rotation of its metric that the + // group lacks (gemmi's find_twin_laws, the same question RfreeFlags and ReindexAmbiguity ask). A + // list of Laue classes cannot answer it - a pseudo-merohedral 2/m or mmm lattice admits a law the + // class names nothing about, and a C-centred orthorhombic lattice with no special metric admits + // none although mmm is not "holohedral" by name. In P1, or without a cell, it cannot be ruled out. + result.merohedral_twinning_possible = + !space_group || !cell || TwinLawPossible(*cell, *space_group); // The two indicators can only move one way under twinning: <|L|> down from 0.500 towards 0.375, // the second moment down from 2.0 towards 1.5. A narrow intensity distribution sitting next to an // <|L|> at or ABOVE its untwinned value therefore has some other cause, and calling it a twin is @@ -356,7 +348,7 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg if (result.merohedral_twinning_possible) result.twinning_suspected = !l_test_contradicts_twin && (l_test_low || moment_low); else { - // In a holohedral class the alternative is not a partial twin anywhere between 0.5 and 0.375 + // With no twin law possible the alternative is not a partial twin anywhere between 0.5 and 0.375 // but a false operator, which puts <|L|> at the perfect-twin 0.375 itself - so the bar sits // nearer that value. Measured on the corpus (normalised L-test): the weakest genuine // holohedral crystal reads 0.436, the over-promoted twin 0.374. The second moment does not @@ -365,6 +357,9 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg // law lies OUTSIDE the Laue class, pulls down (1.66 beside an untwinned <|L|>). result.adopted_operators_suspect = l_test_usable && result.mean_abs_l < 0.42; result.twinning_suspected = result.adopted_operators_suspect; + // Above that bar but still in the twin band: not a twin in this group, which the lattice rules + // out, so not "needs a twin law" - the statistics say instead that the group may be too high. + result.twin_like_without_twin_law = !result.adopted_operators_suspect && l_test_low; } // The twin fraction, from the statistic that carries the verdict. Merged under a false operator @@ -411,8 +406,8 @@ std::string TwinningAnalysisToText(const TwinningAnalysisResult& result) { << " [untwinned 2.00, perfect twin 1.50]\n"; const bool l_test_read = result.l_test_pairs > 0 && !result.l_test_contaminated_by_tncs; if (result.adopted_operators_suspect) - os << " => An adopted operator relates UNEQUAL intensities. The Laue class is holohedral, so no\n" - << " twin law exists within it - yet the intensities read like a twin's, and merging under\n" + os << " => An adopted operator relates UNEQUAL intensities. The lattice admits no twin law for this\n" + << " group - yet the intensities read like a twin's, and merging under\n" << " a false operator produces exactly that distribution, whatever the twin fraction.\n" << (result.laue_class_was_chosen_by_promotion ? " This class was chosen by the space-group search: the promotion is suspect.\n" @@ -440,14 +435,20 @@ std::string TwinningAnalysisToText(const TwinningAnalysisResult& result) { : ".\n"); os << " Statistics flag the presence of twinning, not the twin law; the law itself is a\n" << " question for a dedicated twin-law analysis.\n"; - } else if (!result.merohedral_twinning_possible && result.laue_class_was_chosen_by_promotion) - os << " => No twinning indicated. The Laue class is holohedral, so no twin law exists within it,\n" + } else if (result.twin_like_without_twin_law) + os << " => Twin-like statistics, but the lattice admits no twin law for this group (no rotation of\n" + << " its metric that the group lacks), so they cannot be a twin in this group. Either the\n" + << " symmetry is too high - a twin of a subgroup whose law the group absorbed, or non-\n" + << " crystallographic symmetry taken for real - or the statistics are biased by something\n" + << " other than twinning (non-crystallographic symmetry, overlapping reflections).\n"; + else if (!result.merohedral_twinning_possible && result.laue_class_was_chosen_by_promotion) + os << " => No twinning indicated. The lattice admits no twin law for this group,\n" << " and it was chosen by the space-group search - but its operators read as genuine: a\n" << " merge under a false operator has the perfect-twin distribution, not the untwinned one\n" << " measured here. An operator that is only APPROXIMATE (non-crystallographic symmetry\n" << " along it) also passes; the subgroup statistics reported by the search speak to that.\n"; else if (!result.merohedral_twinning_possible) - os << " => No twinning: the Laue class is holohedral, so no merohedral twin law exists, and no\n" + os << " => No twinning: the lattice admits no twin law for this group, and no\n" << " operator of it averages unequal intensities.\n"; else if (result.l_test_contaminated_by_tncs) os << " => No twinning indicated by the second moment. The L-test, which is normally the\n" @@ -473,7 +474,10 @@ std::string TwinningVerdictLine(const TwinningAnalysisResult& result) { os << std::setprecision(2) << ", twin fraction ~" << result.estimated_twin_fraction << " from the " << (result.twin_fraction_source == TwinFractionSource::LTest ? "L-test" : "second moment"); os << ")"; - } else if (l_test_read) + } else if (result.twin_like_without_twin_law) + os << "TWIN-LIKE, NO TWIN LAW (<|L|> " << result.mean_abs_l + << ": the lattice admits no twin law for this group - the symmetry may be too high)"; + else if (l_test_read) os << "no indication (<|L|> " << result.mean_abs_l << ")"; else os << "no indication (second moment " << std::setprecision(2) << result.second_moment diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index cda213ece..b914d125e 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -22,8 +22,8 @@ // - Second moment /^2 of acentric intensities, normalised per resolution shell so it is // 2.0 for untwinned data and 1.5 for a perfect twin (= 2(1 - a + a^2) for twin fraction a). // -// These flag the *presence* of twinning; they do not by themselves identify the twin law. In a -// holohedral Laue class, where no twin law exists, the same low reading says instead that an operator +// These flag the *presence* of twinning; they do not by themselves identify the twin law. Where the +// lattice admits no twin law for the group, the same low reading says instead that an operator // of the group averages unequal intensities - the group is too high, or a twin law was absorbed into // it (adopted_operators_suspect). @@ -55,11 +55,17 @@ struct TwinningAnalysisResult { // 0.5). Set in every Laue class; what it means depends on the next two fields. bool twinning_suspected = false; - // False when the Laue class is holohedral (4/mmm, 6/mmm, m-3m, rhombohedral -3m): no merohedral - // twin law can exist WITHIN it. + // False when the lattice admits no twin law for the group (TwinLawPossible): no rotation of the + // metric lies outside it, so no merohedral or pseudo-merohedral twin law can exist. bool merohedral_twinning_possible = true; - // twinning_suspected in a holohedral class: an operator of the group averages intensities that + // No twin law is possible, yet <|L|> sits in the twin band (below 0.44) above the bar + // adopted_operators_suspect is read at. Not a twin in this group - the lattice rules that out - so + // a reading that the group may be too high (a subgroup's twin law absorbed, or non-crystallographic + // symmetry taken for real), or that something other than twinning narrowed the distribution. + bool twin_like_without_twin_law = false; + + // twinning_suspected where no twin law is possible: an operator of the group averages intensities that // are not equal, so the group is too high or it has absorbed a twin law. A warning about the // symmetry, not a twin call - and the twin fraction is not measurable on such a merge. bool adopted_operators_suspect = false; @@ -81,10 +87,18 @@ struct TwinningAnalysisResult { bool laue_class_was_chosen_by_promotion = false; }; +// Whether a twin law can exist for `space_group` on this lattice: a rotation of the lattice's own metric +// symmetry - found by Le Page's two-fold search within `max_obliquity_deg`, phenix.xtriage's default of +// 3 degrees, so a pseudo-merohedral lattice counts as its higher metric - that the group lacks. +bool TwinLawPossible(const gemmi::UnitCell& cell, const gemmi::SpaceGroup& space_group, + double max_obliquity_deg = 3.0); + // The space group (when known) is used to drop centric reflections, which follow different // statistics (<|L|>, /^2 = 3 not 2); pass nullptr for P1 data, where none are centric. // `lattice_centring` drops the reflections the lattice centring extinguishes, for a P1 merge of a // centred lattice (whose space_group is nullptr and so says nothing about them). +// `cell` decides whether a twin law is possible at all (TwinLawPossible); without it, or in P1, one is +// assumed possible. // `tncs_vector` is the fractional pseudo-translation of a DETECTED translational pseudo-symmetry on // this same merge (TranslationalNCSVectorForLTest), or nullptr. It changes which partner reflections // the L-test pairs with, so it has to be an argument rather than a flag set on the result afterwards: @@ -93,7 +107,8 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg const gemmi::SpaceGroup* space_group = nullptr, int resolution_shells = 20, const std::array* tncs_vector = nullptr, - char lattice_centring = 'P'); + char lattice_centring = 'P', + const gemmi::UnitCell* cell = nullptr); std::string TwinningAnalysisToText(const TwinningAnalysisResult& result); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 8c3902113..dcd03bc28 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1830,12 +1830,15 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, // NOT_READABLE where <|L|> reads outside the range any twin fraction gives: then neither // twinning nor a false operator is read from it. const bool abnormal = TwinStatisticsAbnormal(result); + // NO_TWIN_LAW: twin-like, but the lattice admits no twin law for the group (see + // TwinningAnalysisResult::twin_like_without_twin_law) - a question about the symmetry. Add(s, KeyEnum("TWINNING_VERDICT", abnormal ? "NOT_READABLE" : tw.adopted_operators_suspect ? "SYMMETRY_SUSPECT" : tw.twinning_suspected ? "INDICATED" + : tw.twin_like_without_twin_law ? "NO_TWIN_LAW" : "NO_INDICATION", - {"NO_INDICATION", "INDICATED", "SYMMETRY_SUSPECT", "NOT_READABLE"})); + {"NO_INDICATION", "INDICATED", "SYMMETRY_SUSPECT", "NO_TWIN_LAW", "NOT_READABLE"})); // The fraction comes from the statistic that carries the verdict, and the key beside it says // which. None is quoted where it cannot be measured: merged under a false operator, or from a // second moment a pseudo-translation inflates. @@ -1862,7 +1865,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, Prose(" Twinning: " + TwinVerdict(result) + ".")); if (!abnormal && tw.adopted_operators_suspect) Add(s, Prose(fmt::format( - " The Laue class is holohedral, so no twin law exists within it - but merging under a\n" + " The lattice admits no twin law for this group - but merging under a\n" " false operator gives exactly this distribution, whatever the twin fraction. {} Either\n" " the crystal has lower symmetry or it is a twin whose law the point group absorbed;\n" " refinement belongs in a subgroup.", @@ -1871,6 +1874,11 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, : "The space group in use is higher than these intensities support."))); else if (!abnormal && tw.twinning_suspected) Add(s, Prose(" Refinement against these data needs a twin law.")); + else if (!abnormal && tw.twin_like_without_twin_law) + Add(s, Prose(" The lattice admits no twin law for this group, so these statistics cannot be a twin in it:\n" + " either the space group is too high - a twin of a subgroup whose law it absorbed, or\n" + " non-crystallographic symmetry taken for real - or something other than twinning\n" + " narrowed the intensity distribution.")); else if (tw.l_test_contaminated_by_tncs) // <|L|> is the number this sentence would otherwise rest on, and on this crystal it is // biased by the pseudo-translation above. @@ -1917,10 +1925,17 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (!abnormal && tw.adopted_operators_suspect) Warn(doc, PathologyCode::TWINNING, fmt::format( - "An adopted symmetry operator relates unequal intensities (<|L|> = {:.3f} in a holohedral " - "Laue class, where a genuine operator reads about 0.5 and a false one 0.375): the " - "space group may be too high, or the crystal is a twin whose law the point group " - "absorbed - refine in a subgroup", + "An adopted symmetry operator relates unequal intensities (<|L|> = {:.3f} where the " + "lattice admits no twin law, so a genuine operator reads about 0.5 and a false one " + "0.375): the space group may be too high, or the crystal is a twin whose law the " + "point group absorbed - refine in a subgroup", + tw.mean_abs_l)); + else if (!abnormal && tw.twin_like_without_twin_law) + Warn(doc, PathologyCode::TWINNING, fmt::format( + "The intensities read twin-like (<|L|> = {:.3f}), but the lattice admits no twin law " + "for this space group, so they cannot be a twin in it: the space group may be too " + "high (a twin of a subgroup whose law it absorbed, or non-crystallographic symmetry " + "taken for real), or something else narrowed the distribution", tw.mean_abs_l)); else if (!abnormal && tw.twinning_suspected) Warn(doc, PathologyCode::TWINNING, fmt::format( diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 65d5b7d04..5a4f04d6f 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -2848,8 +2848,11 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { const char twin_centring = twin_sg ? twin_sg->centring_type() : end_msg.rotation_lattice_type.has_value() ? end_msg.rotation_lattice_type->centering : 'P'; + const std::optional twin_cell = + result.consensus_cell ? std::optional(*result.consensus_cell) : std::nullopt; result.twinning = AnalyzeTwinning(sm.merged, twin_sg, 20, - TranslationalNCSVectorForLTest(result.tncs), twin_centring); + TranslationalNCSVectorForLTest(result.tncs), twin_centring, + twin_cell ? &*twin_cell : nullptr); // Set on the struct the text below is rendered from. AnalyzeTwinning returns a fresh result // with this flag clear, so setting it only afterwards left the statistics text asserting // "No twinning: the Laue class is holohedral" as a fact while the report, rendering the same diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index f370f5be5..6524fb581 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -2150,9 +2150,13 @@ static int RunRugnux(int argc, char **argv) { // Same as the full pipeline: without the reflections the lattice centring extinguishes. The // merged ASU of a centred group still holds keys in the absent class, and they are noise, not // intensities, in both statistics. + const auto twin_uc = experiment.GetUnitCell(); + const std::optional twin_cell = + twin_uc ? std::optional(*twin_uc) : std::nullopt; const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg, 20, TranslationalNCSVectorForLTest(tncs), - twin_sg ? twin_sg->centring_type() : 'P'); + twin_sg ? twin_sg->centring_type() : 'P', + twin_cell ? &*twin_cell : nullptr); logger.Info("{}", TwinningAnalysisToText(twinning)); // Diffraction anisotropy, exactly as the full pipeline computes it (Rugnux.cpp) - no intensity diff --git a/tests/TwinningAnalysisTest.cpp b/tests/TwinningAnalysisTest.cpp index ed51ded0b..165ac5501 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -195,7 +195,8 @@ TEST_CASE("L-test on a P1 merge of a centred lattice", "[twinning]") { // A genuine holohedral crystal: no twin law exists, and its operators read as real. TEST_CASE("No twin is called in a genuine holohedral Laue class", "[twinning]") { const auto merged = MergedIn(Tetragonal("P 4 2 2", 0.0), "P 4 2 2"); - const auto r = AnalyzeTwinning(merged, gemmi::find_spacegroup_by_name("P 4 2 2")); + const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90); + const auto r = AnalyzeTwinning(merged, gemmi::find_spacegroup_by_name("P 4 2 2"), 20, nullptr, 'P', &tetragonal); CHECK_FALSE(r.merohedral_twinning_possible); CHECK(r.mean_abs_l == Catch::Approx(0.5).margin(0.03)); CHECK_FALSE(r.twinning_suspected); @@ -208,8 +209,9 @@ TEST_CASE("No twin is called in a genuine holohedral Laue class", "[twinning]") // law exists". TEST_CASE("A holohedral merge under a false operator is called suspect", "[twinning]") { const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 4 2 2"); + const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90); for (double alpha : {0.0, 0.2, 0.5}) { - const auto r = AnalyzeTwinning(MergedIn(Tetragonal("P 4", alpha), "P 4 2 2"), sg); + const auto r = AnalyzeTwinning(MergedIn(Tetragonal("P 4", alpha), "P 4 2 2"), sg, 20, nullptr, 'P', &tetragonal); CAPTURE(alpha); CHECK(r.mean_abs_l == Catch::Approx(0.375).margin(0.02)); CHECK(r.adopted_operators_suspect); @@ -218,6 +220,53 @@ TEST_CASE("A holohedral merge under a false operator is called suspect", "[twinn } } +// Whether a twin law can exist is read off the lattice, not off a list of Laue classes: a group whose +// lattice has a rotation the group lacks admits one, whatever the class is called. +TEST_CASE("A twin law is possible exactly where the lattice has a rotation the group lacks", "[twinning]") { + const auto possible = [](const char *group, const gemmi::UnitCell &cell) { + return TwinLawPossible(cell, gemmi::get_spacegroup_by_name(group)); + }; + // Merohedral: a tetragonal lattice under a group of lower Laue class than its holohedry. + CHECK(possible("P 4", gemmi::UnitCell(47, 47, 63, 90, 90, 90))); + CHECK_FALSE(possible("P 4 2 2", gemmi::UnitCell(47, 47, 63, 90, 90, 90))); + // Orthorhombic: no law on a general metric - mmm is the lattice's own holohedry - and a law where + // the C-centred cell is pseudo-hexagonal (b = a sqrt(3)). + CHECK_FALSE(possible("C 2 2 21", gemmi::UnitCell(56, 157, 143, 90, 90, 90))); + CHECK(possible("C 2 2 21", gemmi::UnitCell(50, 50 * std::sqrt(3.0), 143, 90, 90, 90))); + // Monoclinic: none on a general metric, one where beta is 90 (pseudo-merohedral). + CHECK_FALSE(possible("P 1 21 1", gemmi::UnitCell(75, 94, 99, 90, 93.1, 90))); + CHECK(possible("P 1 21 1", gemmi::UnitCell(55, 108, 55.3, 90, 90.0, 90))); + // ... and within the obliquity it is read at, a degree off still counts. + CHECK(possible("P 1 21 1", gemmi::UnitCell(54, 142, 84, 90, 91.0, 90))); +} + +// A twin-like <|L|> in a group the lattice admits no twin law for is not a twin in that group: it is +// reported as such, and not as "needs a twin law". Where the lattice does admit one, it is a twin. +TEST_CASE("Twin-like statistics where no twin law is possible", "[twinning]") { + const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 2 2 2"); + const auto merged = MergedIn(Tetragonal("P 2 2 2", 0.12), "P 2 2 2"); + const gemmi::UnitCell orthorhombic(47, 52, 63, 90, 90, 90); + const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90); + + const auto none = AnalyzeTwinning(merged, sg, 20, nullptr, 'P', &orthorhombic); + INFO("<|L|> " << none.mean_abs_l); + REQUIRE(none.mean_abs_l > 0.42); + REQUIRE(none.mean_abs_l < 0.44); + CHECK_FALSE(none.merohedral_twinning_possible); + CHECK_FALSE(none.twinning_suspected); + CHECK_FALSE(none.adopted_operators_suspect); + CHECK(none.twin_like_without_twin_law); + CHECK(TwinningVerdictLine(none).rfind("TWIN-LIKE, NO TWIN LAW", 0) == 0); + + const auto law = AnalyzeTwinning(merged, sg, 20, nullptr, 'P', &tetragonal); + CHECK(law.merohedral_twinning_possible); + CHECK(law.twinning_suspected); + CHECK_FALSE(law.twin_like_without_twin_law); + + // Without a cell nothing can be ruled out. + CHECK(AnalyzeTwinning(merged, sg).merohedral_twinning_possible); +} + namespace { // A P1 merge whose centric reflections are centric: SyntheticMergedReflections draws every // intensity from the acentric distribution, which is all the searches it was written for need, but From 475df9e5cb8d320c9902f1394d22b8861c674370 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 18:59:57 +0200 Subject: [PATCH 3/4] rugnux: write the twin hypothesis beside the answer as _twin_.mtz/.cif Near-perfect twins of a subgroup are undecidable from intensities, and a merge in the higher group cannot be twin-refined. Where the adopted group may have absorbed a twin law - a twin gate refused it and was overruled, or a group above P1 was adopted on a merge that read twin-like before the search - the same observations are also merged in a subgroup and written beside the answer, with the twin law and its fraction in _pdbx_reflns_twin. The subgroup is the one a decision named (the group a refused promotion was adopted over, the indexed lattice's group a metric re-ask promoted from), else the index-2 subgroup the search confirmed as a finalist whose lost two-fold agrees worst with the intensities. It is merged from the same engine, to the answer's resolution, without writing per-frame scales back; it carries the answer's free set (flags read off the answer's merge by the adopted group's asymmetric unit) and is written in the subgroup's reference and standard setting (a C 1 1 2_1 description is reduced to P 1 2_1 1). The fraction is Yeates' H on the pairs the law relates (1/2 - median H, a lower bound). At most one file; other qualifying subgroups are named in the report (TWIN_SUBGROUP_* keys). The answer's files are byte-identical. IndexTwoTwinSubgroups and EstimateTwinLawFraction are new and tested; the point-group ledger now carries each hypothesis's representative group. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/ACKNOWLEDGEMENT.md | 6 + docs/CPU_DATA_ANALYSIS.md | 1 + docs/CPU_DATA_ANALYSIS_DECISIONS.md | 9 +- docs/RUGNUX_REPORT.md | 30 +- docs/RUGNUX_TUTORIAL.md | 4 + .../scale_merge/SearchSpaceGroup.cpp | 2 + image_analysis/scale_merge/SearchSpaceGroup.h | 3 + .../scale_merge/TwinningAnalysis.cpp | 99 ++++++ image_analysis/scale_merge/TwinningAnalysis.h | 26 ++ rugnux/ResultReport.cpp | 28 ++ rugnux/Rugnux.h | 17 + rugnux/RugnuxScaleMerge.cpp | 301 ++++++++++++++++++ tests/TwinningAnalysisTest.cpp | 37 ++- 13 files changed, 557 insertions(+), 6 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index c2de28fed..64064a11d 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -172,6 +172,12 @@ axis preserves the class of a half-integer pseudo-translation, which is what a p That robustness does not extend to a pseudo-translation which is not half-integer, and the partner steps are restricted when one is detected - see the translational-pseudo-symmetry note below. +**The twin fraction of the twin hypothesis** follows Yeates: the intensity a twinned crystal records +is (1−α) I(h) + α I(Th), and the disagreement +H = |I(h) − I(Th)|/(I(h) + I(Th)) of the pairs a law relates is uniform on [0, 1−2α], so the median +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(97)76068-3). + **Translational pseudo-symmetry** is detected from the native Patterson computed from the merged intensities, and the interpretation of an off-origin peak as a pseudo-translation between copies of the contents of the asymmetric unit - together with the modulation it puts on the intensities, which diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index b5b39ccce..b592e9dc5 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -47,6 +47,7 @@ The methods draw on, and in places reimplement, solutions from: - P. Gasparotto, L. Barba, H.-C. Stadler et al., "TORO Indexer: a PyTorch-based indexing algorithm for kilohertz serial crystallography", *J. Appl. Cryst.* **57** (2024), 931-944 (the algorithm of the `ffbidx` fast-feedback indexer, §4). - I. Křivý & B. Gruber, "A unified algorithm for determining the reduced (Niggli) cell", *Acta Cryst.* **A32** (1976), 297-298, and International Tables for Crystallography Vol. A, Table 9.2.5.1 (the Niggli reduction and the lattice-character table of §5.3/§6). - J. E. Padilla & T. O. Yeates, "A statistic for local intensity differences: robustness to anisotropy and pseudo-centering and utility for detecting twinning", *Acta Cryst.* **D59** (2003), 1124-1130 (the L test, §13.2). +- T. O. Yeates, "Detecting and overcoming crystal twinning", *Methods Enzymol.* **276** (1997), 344-358 (the $H$ statistic behind the twin fraction of the twin hypothesis, §13.2). - A. J. C. Wilson, "The probability distribution of X-ray intensities", *Acta Cryst.* **2** (1949), 318-321, and E. R. Howells, D. C. Phillips & D. Rogers, "The probability distribution of X-ray intensities. II. Experimental investigation and the X-ray detection of centres of symmetry", *Acta Cryst.* **3** (1950), 210-214 (the acentric and centric intensity distributions behind the centre-of-symmetry statistics of §13.1 and the Wilson outlier test of §13.3). - W. H. Baur & D. Kassner, "The perils of Cc: comparing the frequencies of falsely assigned space groups with their general population", *Acta Cryst.* **B48** (1992), 356-369 (writing the centrosymmetric group where the absences cannot decide the centre, §13.1). - R. J. Read, P. D. Adams & A. J. McCoy, "Intensity statistics in the presence of translational noncrystallographic symmetry", *Acta Cryst.* **D69** (2013), 176-183 (the native-Patterson detection of translational pseudo-symmetry, and the intensity modulation it produces, which the axial-zone screw-absence test scores against). diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index 122e1f44f..c6812bf84 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -54,13 +54,17 @@ Both search merges also drop the **frames whose fitted per-frame scale came out **A point group's own operators must agree with each other.** A point group is the claim that all of its operators are symmetries of one crystal, and a mixture of real and false operators has a specific shape: the real ones form a subgroup — the true group's intersection with the candidate — and everything outside it is false. The gate therefore reads the operators sorted by their $R$ against the merge, finds every prefix that closes into a subgroup, and takes the largest factor by which every operator outside such a subgroup reads worse than every operator inside it; above 3.2 the candidate is refused. A genuine group whose operators merely spread over a continuum — anisotropic systematic error spreads a cubic crystal's operators twofold — has no such gap. The earlier form, worst-agreeing over best-agreeing operator, divided by an extreme order statistic and refused a genuine $432$ as soon as one of its operators happened to agree unusually well. -**A group that holds every rotation of its lattice is refused when merging under it narrows the L-test like a twin.** Every gate above compares the added operators with a reference, and a pseudo-symmetric structure or a twin defeats them all. For a candidate that holds every rotation its lattice admits no twin law exists, so intensities merged under it that read like a twin's are not a twin. The Padilla–Yeates $L$-test is read twice on the same pairs of the $P1$ merge — as measured, and with every intensity replaced by the mean of its orbit under the candidate — and the candidate is refused where the merged $\langle|L| angle$ reads twin-like ($0.375 \le \langle|L| angle < 0.44$) and the averaging closed a third or more of the gap from the unmerged value to 0.375. Genuine lattice-holohedral groups close 0.01–0.25 of it. A merge reading below 0.375 is not read: nothing a twin or a false operator does reaches that, so something else compresses the intensities. +**A group that holds every rotation of its lattice is refused when merging under it narrows the L-test like a twin.** Every gate above compares the added operators with a reference, and a pseudo-symmetric structure or a twin defeats them all. For a candidate that holds every rotation its lattice admits no twin law exists, so intensities merged under it that read like a twin's are not a twin. The Padilla–Yeates $L$-test is read twice on the same pairs of the $P1$ merge — as measured, and with every intensity replaced by the mean of its orbit under the candidate — and the candidate is refused where the merged $\langle|L| +angle$ reads twin-like ($0.375 \le \langle|L| +angle < 0.44$) and the averaging closed a third or more of the gap from the unmerged value to 0.375. Genuine lattice-holohedral groups close 0.01–0.25 of it. A merge reading below 0.375 is not read: nothing a twin or a false operator does reaches that, so something else compresses the intensities. **The lattice class is re-asked where the metric knows more.** The Bravais class comes from Niggli reduction and a lattice-character lookup (§6), and that lookup can land short of the truth or beside it: near the Niggli type-I/type-II boundary it is decided by the last digits of the refined cell, and a lattice that is nearly but not exactly hexagonal matches the hexagonal character even when no point group that class can host carries the two-folds the data actually have. Two recoveries run, both settled by the intensities. Every rotation the cell metric can host beyond the named class — read off Le Page's two-fold search on the lattice itself, which measures each rotation's obliquity in a primitive basis and owes nothing to the character table — is put to the intensities as a single operator, scored exactly as the search scores its own operators, on the same reflection population and the same $E^2$ normalisation. And where the metric group is larger than the adopted class's holohedry, the merge is reindexed into the metric group's conventional cell and the space-group search is run again there, with every gate live; the reindex is committed only where the search in the new setting confirms a strictly higher point group *and* the centring the new cell describes, so a pseudo-symmetric metric leaves the answer already in hand standing. ### 13.2 Twinning check, and translational pseudo-symmetry -A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500 and 0.333 untwinned, 0.375 and 0.200 for a perfect twin) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (2.0 for untwinned acentric data, 1.5 for a perfect twin) are written to the merged mmCIF as a twinning diagnostic. Both are taken on intensities divided by their resolution-shell mean, over the shells whose $\langle I/\sigma\rangle$ reaches 1, with Wilson outliers rejected; the $L$-test pairs are therefore compared on one scale even where two index steps span a steep fall-off, as in phenix.xtriage and ctruncate, and the selection is by shell, never by the individual reflection's $I/\sigma$, which would cut the weak tail and bias $\langle|L|\rangle$ down. The twin fraction is quoted from the statistic that carries the verdict — the $L$-test unless the call rests on the second moment alone — and a second moment is not turned into a fraction under a detected pseudo-translation, which inflates it. A merohedral twin law exists only where the Laue class is a proper subgroup of the lattice holohedry, so in the high-symmetry holohedral classes ($4/mmm$, $6/mmm$, $m\bar{3}m$, and $\bar{3}m$ on a rhombohedral lattice) no twin is called. The $L$-test is still read there, for a different question: merging $I(h)$ with $I(Th)$ under a false operator $T$ gives $(I(h)+I(Th))/2$ whatever the twin fraction, which has the perfect-twin distribution, so $\langle|L|\rangle$ below 0.42 in a holohedral class is reported as **an adopted operator averaging unequal intensities** — the space group is too high, or a twin law was absorbed into the point group — a warning, never a change to the space group. A genuine operator leaves the untwinned 0.5. Reflections that overlap along a very long axis narrow the distribution the same way, which is one reason it stays a warning. The same numbers measured on the P1 merge the space-group search was given, before any point group was adopted, are reported beside them (`_BEFORE_SEARCH`), with that merge's own pseudo-translation declared and the reflections its lattice centring extinguishes left out. The low-symmetry holohedral classes ($\bar{1}$, $2/m$, $mmm$) also admit no strictly merohedral law, but they stay eligible for the twin call on purpose: *pseudo*-merohedral twinning through an accidentally special metric cannot be ruled out from the symmetry alone, and those are the classes it happens in. +A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500 and 0.333 untwinned, 0.375 and 0.200 for a perfect twin) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (2.0 for untwinned acentric data, 1.5 for a perfect twin) are written to the merged mmCIF as a twinning diagnostic. Both are taken on intensities divided by their resolution-shell mean, over the shells whose $\langle I/\sigma\rangle$ reaches 1, with Wilson outliers rejected; the $L$-test pairs are therefore compared on one scale even where two index steps span a steep fall-off, as in phenix.xtriage and ctruncate, and the selection is by shell, never by the individual reflection's $I/\sigma$, which would cut the weak tail and bias $\langle|L|\rangle$ down. The twin fraction is quoted from the statistic that carries the verdict — the $L$-test unless the call rests on the second moment alone — and a second moment is not turned into a fraction under a detected pseudo-translation, which inflates it. A twin law exists only where the **lattice** has a rotation the group lacks, so whether one is possible is read off the cell, not off the name of the Laue class: the lattice's two-folds are found as for the free set (Le Page, within 3° obliquity — GEMMI's `find_twin_laws`), and where none lies outside the group no twin is called. That covers the holohedral classes ($4/mmm$ on a tetragonal lattice, $6/mmm$, $m\bar{3}m$, $\bar{3}m$ on a rhombohedral one) and equally $mmm$ or $2/m$ on a lattice with no special metric, while a pseudo-merohedral $2/m$ or $mmm$ lattice — a monoclinic $\beta$ near 90°, a C-centred cell with $b \approx a\sqrt{3}$ — admits a law and stays eligible. The $L$-test is still read where no law is possible, for a different question: merging $I(h)$ with $I(Th)$ under a false operator $T$ gives $(I(h)+I(Th))/2$ whatever the twin fraction, which has the perfect-twin distribution, so $\langle|L|\rangle$ below 0.42 there is reported as **an adopted operator averaging unequal intensities** — the space group is too high, or a twin law was absorbed into the point group — a warning, never a change to the space group. A genuine operator leaves the untwinned 0.5. Reflections that overlap along a very long axis narrow the distribution the same way, which is one reason it stays a warning. The same numbers measured on the P1 merge the space-group search was given, before any point group was adopted, are reported beside them (`_BEFORE_SEARCH`), with that merge's own pseudo-translation declared and the reflections its lattice centring extinguishes left out. Between that bar and the twin band's 0.44 the verdict is `NO_TWIN_LAW`: the intensities read twin-like, but the lattice admits no law for the group, so they are not a twin in it — the group may be too high (a subgroup's twin law absorbed, or non-crystallographic symmetry taken for real), or something other than twinning narrowed the distribution. It is reported as such rather than as a twin that "needs a twin law" no refinement program could be given. + +**The twin hypothesis.** Where the adopted group may have absorbed a twin law — a twin gate of the search refused it and was overruled (§13.1), or a group above $P1$ was adopted on a merge that read twin-like before the search — the same observations are also merged in a subgroup and written beside the answer as `_twin_.mtz` / `.cif`, in the subgroup's standard setting, to the answer's resolution and with the answer's free set (each reflection takes the flag of the reflection of the answer's merge its orbit falls in, so twin mates share it). The subgroup is the one a decision named — the group a refused promotion was adopted over, or the indexed lattice's group a metric re-ask promoted from — and otherwise, among the index-2 subgroups of the adopted group that the search confirmed as point-group finalists, the one whose lost two-fold agrees worst with the intensities (lowest operator CC): the most likely twin law. One file is written; the other subgroups that qualified are named in the report. The law is in the mmCIF's `_pdbx_reflns_twin` with its fraction, from Yeates' $H$ on the pairs it relates in the subgroup merge — $H = |I(h)-I(Th)|/(I(h)+I(Th))$ is uniform on $[0, 1-2\alpha]$, so $\alpha = 1/2 - \mathrm{median}\,H$, a lower bound since measurement error only adds to $H$. Nothing the answer is made of depends on it: it is a second merge from the same scaled observations, and the answer's files are byte for byte what they are without it. Refinement programs want twinned data merged in the true lower group, never detwinned, which is why the file is the merge and not a correction of the answer; a near-perfect twin of a subgroup and the higher group predict the same intensities, so which of the two is the crystal is a question for a model. Beside them, on rotation data with `--finalist-ledger`, the run reports **twin-immune zone evidence** for the operators the adopted point group adds over each of its index-2 subgroups, read on a full-resolution merge of the observations in P1 that the option makes for the purpose. Under a twin law $T$, $I_\mathrm{obs}(h)=(1-\alpha)I(h)+\alpha I(Th)$; a reflection whose twin mate is itself up to the subgroup and Friedel is untouched at every $\alpha$, and those are exactly the reflections centric in the group but acentric in the subgroup. They read centric ($\langle|E^2-1|\rangle = 0.968$) if the added operators are real and acentric (0.736) if they are a twin law or a pseudo-symmetry — the one intensity statistic that still separates the two at $\alpha = 0.5$, where every operator statistic reads "real". Each zone is normalised against its own mean in resolution bins (and within the two phase classes of a detected pseudo-translation), read only in shells with $\langle I/\sigma\rangle \ge 5$ because noise inflates every class towards centric, and reported with $n$, a standard error and the centric-over-acentric Wilson log-likelihood ratio in nats (both densities convolved with each reflection's measurement error), beside the acentric control. It is read absolutely, never as the difference to the control: a perfect twin's control (0.541) plus noise inflation would otherwise look like true symmetry. An absolute reading is only as good as the normalisation, though, and the acentric control certifies it: an acentric population reads $-0.130$ nats per reflection when the normalisation is right, and whatever the control reads above that is the normalisation's — anisotropy, a pseudo-translation, a pseudo-centring, noise all inflate every class towards centric alike — so the zone, normalised the same way, carries the same per reflection and the **calibrated** evidence has it taken off (a twinned control reads *below* the expectation, so nothing is taken off a twin). The centric side is read **twinned at the fraction the lattice's other operations show** by their own correlation — the strongest CC of a lattice rotation outside the group, relative to the mean CC of the group's own operators, inverted through $\rho = 2\alpha(1-\alpha)/((1-\alpha)^2+\alpha^2)$ — because a twin by another law reaches a genuine zone exactly as it reaches the control: a genuine 321 crystal twinned by a 622 law read its 2-folds at −206 nats against the untwinned centric density and +209 at that law's fraction of 0.07. The operators' own twin law cannot reach their zone, so the acentric side stays untwinned, and the control's expectation is taken at the same fraction. Measured, a 6/m crystal with a 67 Ų anisotropy read its control at $+0.10$ nats per reflection and the zone of a refused 622 the same, so the zone's $+124$ nats were the normalisation's; calibrated it reads $-115$, and genuine promotions keep $+200$ and above. The calibrated zone is what a refused promotion is decided on (above). It is an in-house method: published practice (Yeates' $H$-test, xtriage) excludes these reflections as uninformative about the twin fraction, and no published test was found that uses them positively; the one precedent here is a trigonal crystal whose twin-immune zone read centric and whose higher group was then confirmed by refinement. @@ -268,3 +272,4 @@ Where a CUDA GPU is present, the refinement's target is evaluated on it: the sam The rotation and translation actually taken are reported (`RIGID_BODY_ROTATION_DEG`, `RIGID_BODY_SHIFT_A`), together with the R-free before it (`R_FREE_BEFORE_RIGID_BODY`), so what the placement bought is visible. The coordinates **as placed** are written as `_model.cif`, and as `_model.pdb` where the PDB format can hold the cell (the fragment-screening tools the file feeds — PanDDA, and the pair dimple produces — take a PDB beside the MTZ), so there is a coordinate file that describes the maps: the input's chains, residues, ligands, waters, B-factors, occupancies and anisotropic $U$s, moved, in the same cell and space group as `.mtz`. Taking the frame from the written reflections rather than from the input model matters — §14.6 may have relabelled them to the model's enantiomorph, which is neither the data's original group nor, necessarily, the model's. It is written **whenever the maps are**, not only where the rigid-body step was committed: the model is re-fractionalized into the data cell and may be relabelled whatever the placement decided, so an unmoved model is still not the input file; and a model the null of §14.5 rejected is scored, placed and mapped like any other, which is precisely the case where the density is worth looking at. + diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index f5f340f83..8634ca78b 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -738,14 +738,38 @@ dropped from the twin verdict **in both directions**: it can no longer indicate longer be read as proof that there is none. The second moment then decides alone. **`TWINNING_VERDICT=`** is `NO_INDICATION`, `INDICATED` (⟨|L|⟩ below 0.44, the phenix.xtriage convention, -or the second moment low, in a Laue class that admits a twin law), `SYMMETRY_SUSPECT` (the same low -reading in a holohedral class, where no twin law exists and a false adopted operator gives the same -distribution) or `NOT_READABLE`. ⟨|L|⟩ has a physical range: 0.5 untwinned, 0.375 a perfect twin, +or the second moment low, where the lattice admits a twin law for the group), `SYMMETRY_SUSPECT` +(⟨|L|⟩ below 0.42 where the lattice admits no twin law - a false adopted operator gives the same +distribution), `NO_TWIN_LAW` (⟨|L|⟩ between 0.42 and 0.44 where the lattice admits no twin law: not a +twin in this group, so the group may be too high, or something other than twinning narrowed the +distribution) or `NOT_READABLE`. Whether a twin law is possible is read off the cell - a rotation of +the lattice's metric, within 3° obliquity, that the group lacks - not off the Laue class, so a +pseudo-merohedral monoclinic or orthorhombic lattice admits one and a general one does not. ⟨|L|⟩ has a physical range: 0.5 untwinned, 0.375 a perfect twin, whatever the law. Where it reads below 0.375, or above 0.55 (on this merge or on the one before the space-group search), no twin fraction explains it - the statistic is distorted, by tNCS, anisotropy, overlapping reflections of a very long cell, or too few reflections - and the verdict is `NOT_READABLE`, with no warning: it says nothing about twinning, nor about the space group. +### The twin hypothesis (`TWIN_SUBGROUP_*`) + +Where the adopted space group may have absorbed a twin law - a twin gate of the space-group search +refused it and was overruled, or a group above P1 was adopted on data that read twin-like before the +search (`L_TEST_MEAN_ABS_L_BEFORE_SEARCH`) - the run also writes the same observations merged in a +subgroup, as `_twin_.mtz` and `.cif` in the subgroup's standard setting, to the same +resolution and with the same free set as `.mtz`, the twin law in the mmCIF's +`_pdbx_reflns_twin`. A near-perfect twin of the subgroup predicts the same intensities as the higher +group, and a refinement program can only twin-refine data merged in the lower group, so the file is +there to refine against with the law where a model prefers it. The answer's files and keys do not +depend on it. At most one is written. + +| key | meaning | +|---|---| +| `TWIN_SUBGROUP_FILE` | The file written, `_twin_.mtz` (and `.cif` beside it) | +| `TWIN_SUBGROUP_SPACE_GROUP` | The subgroup, as written | +| `TWIN_SUBGROUP_LAW` | The twin law on that file's axes, as an hkl triplet | +| `TWIN_SUBGROUP_FRACTION` | The law's fraction from the disagreement of the pairs it relates (1/2 − median H, Yeates' H test): a lower bound | +| `TWIN_SUBGROUP_ALTERNATIVES` | Other confirmed subgroups that qualified; `-S` merges in one | + ## Diffraction anisotropy Section 4 also reports how much the fall-off with resolution depends on **direction**, and whether that is diff --git a/docs/RUGNUX_TUTORIAL.md b/docs/RUGNUX_TUTORIAL.md index 213fe8db5..89eac3cfe 100644 --- a/docs/RUGNUX_TUTORIAL.md +++ b/docs/RUGNUX_TUTORIAL.md @@ -285,6 +285,10 @@ Absorption is corrected only by the empirical surface a rotation run fits and ke per-dataset input layout, and the pair dimple produces — and it is skipped, with a log line, for a cell the PDB format cannot hold. See [Validating against a model](RUGNUX_ADVANCED.md#validating-against-a-model-rugnux---model). +- `_twin_.mtz` / `.cif` — only where the chosen space group may have absorbed a twin law: + the same data merged in a subgroup ``, with the twin law and its fraction, same resolution and + free set as `.mtz`, for refinement with that law. See + [The twin hypothesis](RUGNUX_REPORT.md#the-twin-hypothesis-twin_subgroup_). - `_report.txt` — the **results report**: what the run determined, in a form both a person and a beamline script can read. Always written, next to the files above. See [The results report](RUGNUX_REPORT.md#the-results-report). diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 889e8f635..2a9be66fc 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -2071,6 +2071,8 @@ SearchSpaceGroupResult SearchSpaceGroup( for (const auto& c : pg_cands) { PointGroupLedgerEntry e; e.point_group_hm = c.pg->representative ? c.pg->representative->point_group_hm() : "1"; + if (c.pg->representative) + e.representative = *c.pg->representative; e.order = c.order; e.min_class_cc = c.min_class_cc; e.chi2 = c.chi2; diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 06a28e213..71d0b5ff6 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -724,6 +724,9 @@ struct SearchSpaceGroupOptions { // REPORT-ONLY. Nothing in the search reads these back; the adoption is made exactly as before. struct PointGroupLedgerEntry { std::string point_group_hm; + // The symmorphic group of the hypothesis, in the merge's basis: the name cannot say which axes a + // "2" or a "222" is on, and a caller asking whether a given subgroup is a finalist needs to. + std::optional representative; int order = 0; double min_class_cc = 0.0; // weakest operator CC among the group's own operators double chi2 = std::numeric_limits::quiet_NaN(); // merge chi^2 folding under it diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 01349f9a8..bc88ea5d8 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -19,6 +19,7 @@ #include "gemmi/twin.hpp" #include "../../common/ParallelFor.h" +#include "HKLKey.h" namespace { int64_t PackHKL(int h, int k, int l) { @@ -62,6 +63,104 @@ bool TwinLawPossible(const gemmi::UnitCell& cell, const gemmi::SpaceGroup& space return !gemmi::find_twin_laws(cell, &space_group, max_obliquity_deg, /*all_ops=*/false).empty(); } +std::vector IndexTwoTwinSubgroups(const gemmi::SpaceGroup& group) { + const gemmi::GroupOps gops = group.operations(); + const size_t order = gops.sym_ops.size(); + std::vector rots; + for (const auto& op : gops.sym_ops) + rots.push_back(op.rot); + const auto index_of = [&](const gemmi::Op::Rot& r) { + return static_cast(std::find(rots.begin(), rots.end(), r) - rots.begin()); + }; + const auto product = [](const gemmi::Op::Rot& a, const gemmi::Op::Rot& b) { + return gemmi::Op{a, {0, 0, 0}, 'x'}.combine(gemmi::Op{b, {0, 0, 0}, 'x'}).rot; + }; + // Every index-2 subgroup of a crystallographic point group is generated by at most two of its + // operations, so the closures of all pairs find them all. + std::vector> seen; + std::vector out; + for (size_t i = 0; i < order; ++i) + for (size_t j = i; j < order; ++j) { + std::vector in(order, false); + in[0] = in[i] = in[j] = true; // sym_ops[0] is the identity + for (bool grew = true; grew;) { + grew = false; + for (size_t a = 0; a < order; ++a) + for (size_t b = 0; b < order; ++b) + if (in[a] && in[b]) { + const size_t c = index_of(product(rots[a], rots[b])); + if (c < order && !in[c]) { in[c] = true; grew = true; } + } + } + if (2 * static_cast(std::count(in.begin(), in.end(), true)) != order + || std::find(seen.begin(), seen.end(), in) != seen.end()) + continue; + seen.push_back(in); + gemmi::GroupOps sub; + sub.cen_ops = gops.cen_ops; + std::vector laws; + for (size_t a = 0; a < order; ++a) { + if (in[a]) { + sub.sym_ops.push_back(gops.sym_ops[a]); + continue; + } + const gemmi::Op r{rots[a], {0, 0, 0}, 'x'}; + if (r.det_rot() > 0 && r.combine(r).rot == gemmi::Op::identity().rot) + laws.push_back(r); + } + if (laws.empty()) + continue; + // The kept screws can sit off the origin GEMMI's table puts them at (the 2_1 axes of + // P2_12_12_1 do not meet), so the subgroup is looked up at every origin shift; a shift moves + // no reflection and no absence. + const gemmi::SpaceGroup* found = nullptr; + for (int n = 0; n < 12 * 12 * 12 && found == nullptr; ++n) { + gemmi::GroupOps shifted = sub; + gemmi::Op shift = gemmi::Op::identity(); + shift.tran = {n / 144 * gemmi::Op::DEN / 12, n / 12 % 12 * gemmi::Op::DEN / 12, + n % 12 * gemmi::Op::DEN / 12}; + shifted.change_basis_forward(shift); + found = gemmi::find_spacegroup_by_ops(shifted); + } + if (found != nullptr) + out.push_back({*found, laws}); + } + return out; +} + +// Following Yeates (1997) Methods Enzymol. 276, 344-358 +TwinLawFraction EstimateTwinLawFraction(const std::vector& merged, + const gemmi::SpaceGroup& space_group, const gemmi::Op& law, + bool merge_friedel) { + TwinLawFraction result; + const HKLKeyGenerator key_of(merge_friedel, space_group); + const gemmi::GroupOps gops = space_group.operations(); + std::unordered_map intensity; + for (const auto& r : merged) + if (UsableIntensity(r)) + intensity.emplace(key_of(r).pack(), r.I); + std::vector h; + for (const auto& r : merged) { + if (!UsableIntensity(r) || gops.is_reflection_centric({{r.h, r.k, r.l}})) + continue; + const auto mate = law.apply_to_hkl({{r.h, r.k, r.l}}); + const uint64_t k1 = key_of(r).pack(), k2 = key_of(mate[0], mate[1], mate[2]).pack(); + if (k2 <= k1) + continue; // each pair once, and not a reflection with itself + const auto it = intensity.find(k2); + if (it == intensity.end() || r.I + it->second <= 0.0) + continue; + h.push_back(std::fabs(r.I - it->second) / (r.I + it->second)); + } + result.pairs = static_cast(h.size()); + if (h.empty()) + return result; + const size_t mid = h.size() / 2; + std::nth_element(h.begin(), h.begin() + mid, h.end()); + result.fraction = std::max(0.0, 0.5 - h[mid]); + return result; +} + TwinningAnalysisResult AnalyzeTwinning(const std::vector& merged, const gemmi::SpaceGroup* space_group, int resolution_shells, const std::array* tncs_vector, char lattice_centring, diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index b914d125e..66211c0fa 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -4,6 +4,7 @@ #pragma once #include +#include #include #include @@ -112,6 +113,31 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg std::string TwinningAnalysisToText(const TwinningAnalysisResult& result); +// An index-2 subgroup of `group` that keeps the translations of the operations it keeps (a +// translationengleiche subgroup, in the group's own setting), with the proper two-folds of `group` +// outside it: each is a twin law that merge in the subgroup would carry, were `group` a twin's +// supergroup, and any one of them describes the twin (they differ by an operation of the subgroup). +// Only subgroups whose coset holds a two-fold and that GEMMI's table names in this setting (at some origin) are listed. +// The laws are rotations (zero translation), acting on Miller indices through apply_to_hkl. +struct TwinSubgroup { + gemmi::SpaceGroup subgroup; + std::vector laws; +}; +std::vector IndexTwoTwinSubgroups(const gemmi::SpaceGroup& group); + +// 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 +// fraction a, so a = 1/2 - median H. A lower bound: measurement error only adds to H. NaN with no +// pairs. +struct TwinLawFraction { + double fraction = std::numeric_limits::quiet_NaN(); + int pairs = 0; +}; +TwinLawFraction EstimateTwinLawFraction(const std::vector& merged, + const gemmi::SpaceGroup& space_group, const gemmi::Op& law, + bool merge_friedel = true); + // The verdict in one line - "INDICATED (...)", "SYMMETRY SUSPECT (...)" or "no indication (...)" - // so the report's summary row, its prose and the viewer all say the same thing. std::string TwinningVerdictLine(const TwinningAnalysisResult& result); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index dcd03bc28..497d3150b 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -1906,6 +1906,34 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, " against these data in this space group with that law.", result.twinning.twin_law_operator, result.twinning.twin_law_fraction))); } + if (result.twin_subgroup_merge) { + const auto &t = *result.twin_subgroup_merge; + Add(s, KeyText("TWIN_SUBGROUP_FILE", t.file_prefix + ".mtz")); + Add(s, KeyText("TWIN_SUBGROUP_SPACE_GROUP", t.space_group)); + Add(s, KeyText("TWIN_SUBGROUP_LAW", t.twin_law)); + Add(s, KeyReal("TWIN_SUBGROUP_FRACTION", t.twin_fraction, "{:.3f}")); + Add(s, KeyInt("TWIN_SUBGROUP_FRACTION_PAIRS", t.twin_fraction_pairs, true)); + if (!t.other_subgroups.empty()) { + std::string others; + for (const auto &o : t.other_subgroups) + others += (others.empty() ? "" : " | ") + o; + Add(s, KeyText("TWIN_SUBGROUP_ALTERNATIVES", others)); + } + Add(s, Prose(fmt::format( + " Twin hypothesis: {} was adopted, but {}. A twin of {} whose law {} absorbs predicts\n" + " these intensities too, and a merge in {} cannot be refined against that twin. So the\n" + " same observations are also written merged in {} - same resolution, same free set -\n" + " as {}.mtz / .cif, with the twin law {} at a fraction of at least {:.2f} (from the\n" + " disagreement of the reflection pairs it relates, a lower bound; the mmCIF carries it\n" + " as _pdbx_reflns_twin). Refine against it with that twin law where the model prefers it.\n" + " Nothing in this report or in the answer's files depends on it.{}", + t.adopted_space_group, t.reason, t.space_group, t.adopted_space_group, t.space_group, + t.space_group, t.file_prefix, t.twin_law, t.twin_fraction, + t.other_subgroups.empty() ? std::string() + : "\n Other subgroups the search confirmed, each a twin hypothesis of\n" + " its own (re-run with -S to merge in one): see\n" + " TWIN_SUBGROUP_ALTERNATIVES."))); + } for (size_t i = 0; i < result.twin_law_evidence.size(); ++i) Add(s, KeyText(fmt::format("TWIN_LAW_{}", i + 1).c_str(), result.twin_law_evidence[i])); if (result.twin_immune_zones && !result.twin_immune_zones->zones.empty()) { diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 4b5c3dd46..c6eb42c66 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -477,6 +477,23 @@ struct ProcessResult { // merge (AnalyzeTwinImmuneZones). Report-only. Unset without --finalist-ledger. std::optional twin_immune_zones; + // The twin hypothesis written beside the answer: where the adopted group may have absorbed a twin + // law (a twin gate overruled, or a group adopted on data that read twin-like before the search), the + // same observations merged in a subgroup, to the same resolution and with the same free set, as + // _twin_.mtz / .cif in the subgroup's standard setting, carrying the law. The answer's + // own files are not touched by it. Unset where none was written. + struct TwinSubgroupMerge { + std::string reason; // why the hypothesis was written + std::string adopted_space_group; // the group of the answer, as written + std::string space_group; // the subgroup, as written (standard setting) + std::string file_prefix; // _twin_ + std::string twin_law; // hkl triplet, on the subgroup file's axes + double twin_fraction = NAN; // EstimateTwinLawFraction: a lower bound + int twin_fraction_pairs = 0; + std::vector other_subgroups; // other confirmed subgroups that qualified + }; + std::optional twin_subgroup_merge; + // What --model found, when one was given. Set even when the validation could not run - its // failure_reason then says why - so the report can report the failure rather than leave a run // that was asked to validate looking like one that was not. Unset when no model was asked for. diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 5a4f04d6f..12641df51 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -214,6 +214,108 @@ namespace { } return best; } + + // The subgroup a twin-law merge is made in (_twin_), and the law it carries, where the + // adopted group may have absorbed a twin law. Candidates are the index-2 subgroups of `group` + // (IndexTwoTwinSubgroups); where a decision of the search named one - the group a refused promotion + // was adopted over, or the indexed lattice's group a metric re-ask promoted from - only that one, + // and otherwise those the search confirmed as point-group finalists (the ledger's eligible rows). + // Among them, the one whose law agrees worst with the intensities (lowest Stage A CC): the + // two-fold most likely to be the twin law. The rest are named for the report. + struct TwinSubgroupChoice { + gemmi::SpaceGroup subgroup; + gemmi::Op law; + std::vector others; + }; + + std::optional ChooseTwinSubgroup(const gemmi::SpaceGroup &group, + const SearchSpaceGroupResult &search, + const std::optional &named) { + 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 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(); + }; + std::optional best; + double best_cc = 0.0; + std::vector names; + for (const auto &t : IndexTwoTwinSubgroups(group)) { + if (t.subgroup.operations().sym_ops.size() < 2) + continue; // P1: a twin of it is not what a promotion absorbs + const auto rots = rotations(t.subgroup); + const bool qualifies = named + ? rotations(*named) == rots + : std::any_of(search.point_group_ledger.begin(), search.point_group_ledger.end(), + [&](const PointGroupLedgerEntry &e) { + return e.eligible && e.representative + && rotations(*e.representative) == rots; + }); + if (!qualifies) + continue; + names.push_back(t.subgroup.xhm()); + const gemmi::Op *law = &t.laws.front(); + for (const auto &l : t.laws) + if (law_cc(l) < law_cc(*law)) + law = &l; + if (!best || law_cc(*law) < best_cc) { + best = TwinSubgroupChoice{t.subgroup, *law, {}}; + best_cc = law_cc(*law); + } + } + if (best) + for (const auto &n : names) + if (n != best->subgroup.xhm()) + best->others.push_back(n); + return best; + } + + // 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. + bool MoveTwinSubgroupChoice(TwinSubgroupChoice &c, const gemmi::Op &cob) { + const gemmi::SpaceGroup *moved = SpaceGroupInBasis(c.subgroup, cob); + if (moved == nullptr) + return false; + const gemmi::Op m = HklOperator(cob); + c.subgroup = *moved; + c.law = m.inverse().combine(c.law).combine(m); // h' = m h: back by m, the law, forward by m + return true; + } + + // `group` carried into the setting a Miller-index reindex q (h' = q h) takes the merge to, with the + // centring that setting describes - the move the metric re-ask makes with the lattice. Empty where + // no setting in GEMMI's table names the result. + std::optional GroupInReindexedSetting(const gemmi::SpaceGroup &group, + const gemmi::Mat33 &q, char centring) { + try { + gemmi::Op hkl = gemmi::Op::identity(); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + hkl.rot[i][j] = static_cast(std::lround(q[j][i])) * gemmi::Op::DEN; // apply_to_hkl uses rot^T + gemmi::Op cob = hkl.inverse(); + for (int i = 0; i < 12; ++i) + for (int j = 0; j < 12; ++j) + for (int k = 0; k < 12; ++k) { + cob.tran = {i * gemmi::Op::DEN / 12, j * gemmi::Op::DEN / 12, k * gemmi::Op::DEN / 12}; + gemmi::GroupOps gops = group.operations(); + gops.change_basis_forward(cob); + gops.cen_ops = gemmi::centring_vectors(centring); + if (const gemmi::SpaceGroup *found = gemmi::find_spacegroup_by_ops(gops)) + return *found; + } + } catch (const std::exception &) { + } + return std::nullopt; + } } // Proper rotations of a crystal system's holohedry - what a lattice of that class carries, and the @@ -1210,6 +1312,16 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // above triclinic - i.e. a promotion was made, which is what makes a later "the Laue class is // holohedral so there is no twin law" conclusion circular. bool promoted_point_group = false; + // The twin-law merge (_twin_): where the adopted group may have absorbed a twin law, + // the subgroup a second merge is made in and the law it carries, in the scaler's frame, and why. + // The decisions below that name a subgroup leave it here: the group a refused promotion was + // adopted over (a twin gate overruled), and the indexed lattice's group a metric re-ask promoted + // from. Report and an extra file only - nothing the run writes as its answer reads any of it. + std::optional twin_overruled_subgroup, twin_metric_subgroup; + std::string twin_overruled_reason; + std::optional twin_choice; + std::string twin_choice_reason; + int twin_choice_group_number = 0; // the adopted group's, when the choice was made TimingMark search_mark(logger, "space-group search"); if (search_space_group) { SearchSpaceGroupOptions sg_opts; @@ -1789,6 +1901,20 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { "group already found" : " - keeping the group already found; " + s2.refused_reason))); if (higher) { + // The group of the lattice in hand, carried into the new setting: the subgroup of + // the twin-law merge, should this promotion stand on twin-like data or over a + // refusal by a twin gate (see twin_choice). + if (sg_search.best_space_group) { + gemmi::Mat33 q = reindex; + if (commit_reindex) + q = reindex.multiply(commit_reindex->inverse()); + twin_metric_subgroup = GroupInReindexedSetting(*sg_search.best_space_group, q, + cand.centering); + if (sg_search.refused_by_twin_gate) { + twin_overruled_subgroup = twin_metric_subgroup; + twin_overruled_reason = sg_search.refused_reason; + } + } sg_search = s2; // Leave the candidate it displaced among the alternatives, so the reported // "or ..." list stays the set of groups the data cannot separate. @@ -2476,6 +2602,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { sg_search.refused_point_group_hm, sg_search.refused_reason, hi_sg->xhm()); logger.Info("{}", msg); stats_text << " " << msg << "\n\n"; + if (sg_search.refused_by_twin_gate) { + twin_overruled_subgroup = sg; + twin_overruled_reason = sg_search.refused_reason; + } sg_search = hi_search; if (hi_search.best_space_group->number != hi_sg->number) { sg_search.alternatives.push_back(*hi_search.best_space_group); @@ -2537,6 +2667,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { measured, hi_sg->xhm()); logger.Info("{}", msg); stats_text << " " << msg << "\n\n"; + if (sg_search.refused_by_twin_gate) { + twin_overruled_subgroup = sg; + twin_overruled_reason = sg_search.refused_reason; + } sg_search = hi_search; // Leave the candidate it displaced among the alternatives, so the // reported "or ..." list stays the set of groups the data cannot @@ -2577,6 +2711,28 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { result.space_group.value_or(gemmi::get_spacegroup_p1())); } result.space_group_search = sg_search; + + // Was a twin law possibly absorbed? Where a twin gate refused the adopted group and was + // overruled, or where a group above P1 was adopted on data that read twin-like before the + // search (the pre-promotion L-test), the merge in the subgroup is written beside the answer + // too, with the law - so a refinement program can take the twin hypothesis without a re-run. + // The answer itself is not touched by any of this. + const bool pre_twin_like = result.pre_promotion_twinning && result.pre_promotion_twinning->twinning_suspected; + if (!geometry_prepass && result.space_group && sg_search.point_group_order > 1 + && (twin_overruled_subgroup || pre_twin_like)) { + twin_choice = ChooseTwinSubgroup(*result.space_group, sg_search, + twin_overruled_subgroup ? twin_overruled_subgroup + : twin_metric_subgroup); + twin_choice_group_number = result.space_group->number; + twin_choice_reason = twin_overruled_subgroup + ? fmt::format("a twin gate refused it ({}) and was overruled", twin_overruled_reason) + : fmt::format("the merge read twin-like before the space-group search ({})", + TwinningVerdictLine(*result.pre_promotion_twinning)); + if (twin_choice) + logger.Info("Twin hypothesis: {} was adopted, but {} - the merge in {} is written too, " + "with twin law {}", result.space_group->xhm(), twin_choice_reason, + twin_choice->subgroup.xhm(), twin_choice->law.as_hkl().triplet()); + } } else { // A space group was fixed by the user; surface it so the viewer/CLI can still show it. result.space_group = experiment_.GetGemmiSpaceGroup(); @@ -2701,6 +2857,10 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { if (!(choice.cob == gemmi::Op::identity())) reseat_integrated(choice.cob); relabel_group(choice.cob, *moved); + // The engine is about to be remade in the reference's frame, and the twin-law merge is + // made from it: its subgroup moves with the data. + if (twin_choice && !MoveTwinSubgroupChoice(*twin_choice, choice.cob)) + twin_choice.reset(); determined_to_written = determined_to_written.combine(hkl); result.setting_operator = IndexTriplet(determined_to_written); rsm.emplace(experiment_, indexer->GetIntegrationOutcome(), result.consensus_cell, @@ -2776,6 +2936,24 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // of its own rather than the experiment's, and writes no per-frame scale back onto the outcomes; // and nothing up to take_ledger_p1 calls the engine, changes the outcomes it was ingested from, // or changes the experiment. Its log lines are held and printed where it is taken. + // The twin-law merge (see twin_choice), from the same engine as the answer's, to the resolution the + // answer was cut at. Like the ledger's below, it merges the engine's own copy of the observations in + // a group of its own and writes no per-frame scale back, so nothing the answer is made of sees it. + std::optional twin_merge; + if (twin_choice && rsm && result.consensus_cell && write_files && config_.write_merged + && !geometry_prepass && !superseded) { + TimingMark mark(logger, "twin-law merge"); + phase("Merging in " + twin_choice->subgroup.short_name() + " (twin hypothesis)"); + rsm->SetWriteBackPerFrameScale(false); + rsm->SetExportScaledFulls(false); + rsm->SetFrenchWilson(false); + const auto manual_d_min = experiment_.GetScalingSettings().GetHighResolutionLimit_A(); + rsm->SetDMinLimit(sm.statistics.overall.d_min); + twin_merge = rsm->Run(/*for_search=*/false, /*full_stats=*/true, + /*measure_cc_before_corrections=*/false, &twin_choice->subgroup); + rsm->SetDMinLimit(manual_d_min); + rsm->SetWriteBackPerFrameScale(true); + } std::optional ledger_p1; Logger ledger_log = Logger::Buffered(); std::future ledger_ahead; @@ -3271,6 +3449,109 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // measurements, so the statistics computed above are unaffected; what changes is that the file, // the R-factors and the maps then all describe one indexing instead of two. std::future placed_model_file; + + // The twin-law merge made ready to be written: on the axes the answer is written on (every + // relabelling since it was merged), with the answer's own free set - read off the answer's merge + // by the adopted group's asymmetric unit, which both twin mates of a pair fall in - then in the + // subgroup's standard setting, with its French-Wilson amplitudes and the law's fraction. Made once + // the model, which can relabel the answer, is done with it; and before the model is asked about + // the hypothesis itself (below). + struct TwinFile { + std::vector merged; + UnitCell cell; + gemmi::SpaceGroup space_group; + gemmi::Op law; + TwinLawFraction fraction; + }; + std::optional twin_file; + const auto prepare_twin_file = [&] { + if (twin_file || !twin_merge || !twin_choice) + return; + const bool friedel = experiment_.GetScalingSettings().GetMergeFriedel(); + const gemmi::SpaceGroup adopted = experiment_.GetSpaceGroupOrP1(); + TwinSubgroupChoice c = *twin_choice; + std::vector merged = twin_merge->merged; + if (!(merge_to_written == gemmi::Op::identity())) { + if (!MoveTwinSubgroupChoice(c, merge_to_written.inverse())) { + logger.Warning("Twin hypothesis: {} has no setting on the axes the answer is written on; " + "its merge is not written", twin_choice->subgroup.xhm()); + return; + } + merged = ReindexMergedIntoAsu(merged, merge_to_written, c.subgroup, friedel); + } + // A subgroup of the group written - which the model can relabel to its enantiomorph, a change + // of label only, so the subgroup's label follows it. + const auto rotations_in = [](const gemmi::SpaceGroup &sub, const gemmi::SpaceGroup &group) { + const gemmi::GroupOps g = group.operations(); + for (const auto &op : sub.operations().sym_ops) + if (std::none_of(g.sym_ops.begin(), g.sym_ops.end(), + [&](const gemmi::Op &o) { return o.rot == op.rot; })) + return false; + return true; + }; + if (adopted.number != twin_choice_group_number && c.subgroup.is_enantiomorphic()) { + gemmi::GroupOps e = c.subgroup.operations(); + e.change_basis_forward(c.subgroup.change_of_hand_op()); + if (const gemmi::SpaceGroup *other = gemmi::find_spacegroup_by_ops(e)) + c.subgroup = *other; + } + if (!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; + } + AssignRfreeFlags(merged, c.subgroup, experiment_.GetScalingSettings().GetRfreeFraction(), 500, + result.consensus_cell, config_.nthreads); + ApplyReferenceFreeFlags(merged, adopted, sm.merged); + const TwinLawFraction fraction = EstimateTwinLawFraction(merged, c.subgroup, c.law, friedel); + UnitCell cell = *result.consensus_cell; + // The subgroup's reference setting first: a subgroup of a centred group can be a primitive + // group written on a multiple of its cell - C 1 1 2_1 is P 1 2_1 1 on twice the volume - which + // no unimodular relabelling reaches. The reflections the centring extinguished have no index on + // the primitive cell and are dropped; the rest move exactly. + if (!(c.subgroup.basisop() == gemmi::Op::identity())) { + const gemmi::Op to_reference = c.subgroup.basisop().inverse(); + const gemmi::Op m = HklOperator(to_reference); + const gemmi::SpaceGroup before = c.subgroup; + if (MoveTwinSubgroupChoice(c, to_reference)) { + std::erase_if(merged, [&](const MergedReflection &r) { + const auto h = m.apply_to_hkl_without_division({{r.h, r.k, r.l}}); + return h[0] % gemmi::Op::DEN != 0 || h[1] % gemmi::Op::DEN != 0 || h[2] % gemmi::Op::DEN != 0; + }); + merged = ReindexMergedIntoAsu(merged, m, c.subgroup, friedel); + cell = CellInBasis(cell, to_reference); + } else { + c.subgroup = before; + } + } + const gemmi::Op standard = ChooseOutputSetting(cell, c.subgroup); + if (!(standard == gemmi::Op::identity())) { + const gemmi::SpaceGroup before = c.subgroup; + if (MoveTwinSubgroupChoice(c, standard)) { + merged = ReindexMergedIntoAsu(merged, HklOperator(standard), c.subgroup, friedel); + cell = CellInBasis(cell, standard); + } else { + c.subgroup = before; + } + } + FrenchWilsonOptions fw_opts; + fw_opts.num_threads = config_.nthreads; + ApplyFrenchWilson(merged, c.subgroup, fw_opts); + twin_file = TwinFile{std::move(merged), cell, c.subgroup, c.law, fraction}; + + auto &t = result.twin_subgroup_merge.emplace(); + t.reason = twin_choice_reason; + t.adopted_space_group = adopted.xhm(); + t.space_group = c.subgroup.xhm(); + t.file_prefix = config_.output_prefix + "_twin_" + c.subgroup.short_name(); + t.twin_law = c.law.as_hkl().triplet(); + t.twin_fraction = fraction.fraction; + t.twin_fraction_pairs = fraction.pairs; + t.other_subgroups = c.others; + 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_prefix); + }; if (result.consensus_cell && write_files && !config_.model_path.empty() && !geometry_prepass && !superseded) { phase("Validating against model"); @@ -3464,6 +3745,7 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { WritePlacedModel(*placed, cell, sg, prefix, logger); }); } + prepare_twin_file(); // The unmerged MTZ (written with the other output files, see there) is built from here on, // beside the reflection files: the indexing and the group it is written in are final once the @@ -3508,6 +3790,25 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, em_report, result.twinning, config_.provenance, config_.output_prefix, static_cast(std::max(1, config_.nthreads)), &sm.scaled_fulls); + // The twin hypothesis beside it (see prepare_twin_file): the reflection files only, with the + // law and its fraction in the mmCIF's _pdbx_reflns_twin. + if (twin_file) { + DiffractionExperiment twin_x = experiment_; + twin_x.SetSpaceGroup(twin_file->space_group); + TwinningAnalysisResult law; + law.twin_law_operator = twin_file->law.as_hkl().triplet(); + law.twin_law_fraction = std::isfinite(twin_file->fraction.fraction) ? twin_file->fraction.fraction : 0.5; + const ErrorModelReport twin_em{ + twin_merge->isa > 0 && twin_merge->isa_resolved ? fmt::format("{:.2f}", twin_merge->isa) : "?", + std::string(), twin_merge->error_model_a > 0 ? fmt::format("{:.3f}", twin_merge->error_model_a) + : std::string(), + twin_merge->error_model_b > 0 ? fmt::format("{:.4e}", twin_merge->error_model_b) + : std::string()}; + const std::string prefix = result.twin_subgroup_merge->file_prefix; + WriteMtzReflections(twin_file->merged, twin_file->cell, twin_x, config_.provenance, prefix + ".mtz"); + WriteMmcifReflections(twin_file->merged, twin_file->cell, twin_x, twin_merge->statistics, twin_em, + law, config_.provenance, prefix + ".cif", static_cast(std::max(1, config_.nthreads))); + } writing_mark.Stop(); // FINALIST LEDGER (--finalist-ledger), report-only and deliberately behind a flag: these diff --git a/tests/TwinningAnalysisTest.cpp b/tests/TwinningAnalysisTest.cpp index 165ac5501..92949a622 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -3,7 +3,9 @@ #include +#include #include +#include #include #include #include @@ -244,7 +246,7 @@ TEST_CASE("A twin law is possible exactly where the lattice has a rotation the g // reported as such, and not as "needs a twin law". Where the lattice does admit one, it is a twin. TEST_CASE("Twin-like statistics where no twin law is possible", "[twinning]") { const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 2 2 2"); - const auto merged = MergedIn(Tetragonal("P 2 2 2", 0.12), "P 2 2 2"); + const auto merged = MergedIn(Tetragonal("P 2 2 2", 0.15), "P 2 2 2"); const gemmi::UnitCell orthorhombic(47, 52, 63, 90, 90, 90); const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90); @@ -448,3 +450,36 @@ TEST_CASE("L-test under a merge tells a false operator from a real one", "[twinn CHECK(perfect_twin.mean_abs_l == Catch::Approx(0.375).margin(0.02)); CHECK(perfect_twin.mean_abs_l_merged == Catch::Approx(perfect_twin.mean_abs_l).margin(0.01)); } + +// The subgroups a twin-law merge can be made in: index 2, translations kept, and a two-fold of the +// group outside them to be the law. A four-fold over its two-fold offers no two-fold, so no law. +TEST_CASE("Index-2 subgroups with a two-fold twin law", "[twinning]") { + const auto names = [](const char *group) { + std::vector out; + for (const auto &t : IndexTwoTwinSubgroups(gemmi::get_spacegroup_by_name(group))) { + CHECK_FALSE(t.laws.empty()); + out.push_back(t.subgroup.xhm()); + } + std::sort(out.begin(), out.end()); + return out; + }; + CHECK(names("P 2 2 2") == std::vector{"P 1 1 2", "P 1 2 1", "P 2 1 1"}); + CHECK(names("P 21 21 21") == std::vector{"P 1 1 21", "P 1 21 1", "P 21 1 1"}); + CHECK(names("P 4").empty()); + const auto p422 = names("P 4 2 2"); + CHECK(std::find(p422.begin(), p422.end(), "P 4") != p422.end()); + CHECK(std::find(p422.begin(), p422.end(), "P 2 2 2") != p422.end()); + CHECK(names("C 2 2 21").size() == 3); +} + +// 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"); + const gemmi::Op law = jfjoch_test::TwinLaw(p4, gemmi::get_spacegroup_by_name("P 4 2 2")); + for (double alpha : {0.0, 0.2, 0.4}) { + CAPTURE(alpha); + const auto f = EstimateTwinLawFraction(MergedIn(Tetragonal("P 4", alpha), "P 4"), p4, law); + CHECK(f.pairs > 1000); + CHECK(f.fraction == Catch::Approx(alpha).margin(0.03)); + } +} From 8442d55f6f95dee2b9b0e5a147b8842ba8588347 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 9 Oct 2026 19:02:57 +0200 Subject: [PATCH 4/4] rugnux --model: score the twin hypothesis too, and take the maps from the merge the model prefers Where a twin hypothesis (_twin_) was written, the model is also validated against that merge: as merged, and as twinned by its law - |F_twin(h)|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 from the placed, scaled model, one overall scale refitted per a, a scanned 0..0.5 on the working set and R-free read at it (ValidateAgainstModel's new twin_law argument; GPU structure factors as before). The maps and the placed model are then written from whichever merge has the lower R-free; the report gives MAPS_FROM and the subgroup's plain and twinned R-work/R-free and fraction. The subgroup file follows its own validation's indexing, as the answer follows its own. The model chooses the maps only: the space group, .mtz and R_FREE stay the answer's. Documented in RUGNUX_REPORT.md, RUGNUX_ADVANCED.md and CPU_DATA_ANALYSIS_DECISIONS.md (14.9); Yeates (1997) credited. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01SVmAWnzCmRKAXVUCdc4iNi --- docs/ACKNOWLEDGEMENT.md | 4 +- docs/CPU_DATA_ANALYSIS.md | 2 +- docs/CPU_DATA_ANALYSIS_DECISIONS.md | 8 ++- docs/RUGNUX_ADVANCED.md | 10 ++- docs/RUGNUX_REPORT.md | 8 +++ docs/RUGNUX_TUTORIAL.md | 3 +- .../structure_refinement/ModelValidation.cpp | 64 ++++++++++++++++++- .../structure_refinement/ModelValidation.h | 14 +++- rugnux/ResultReport.cpp | 27 ++++++++ rugnux/Rugnux.h | 4 ++ rugnux/RugnuxScaleMerge.cpp | 56 +++++++++++++++- 11 files changed, 189 insertions(+), 11 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 64064a11d..91ca925f4 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -172,8 +172,8 @@ axis preserves the class of a half-integer pseudo-translation, which is what a p That robustness does not extend to a pseudo-translation which is not half-integer, and the partner steps are restricted when one is detected - see the translational-pseudo-symmetry note below. -**The twin fraction of the twin hypothesis** follows Yeates: the intensity a twinned crystal records -is (1−α) I(h) + α I(Th), and the disagreement +**The twin fraction of the twin hypothesis, and the twinned model it is scored with**, follow Yeates: +the intensity a twinned crystal records is (1−α) I(h) + α I(Th), and the disagreement H = |I(h) − I(Th)|/(I(h) + I(Th)) of the pairs a law relates is uniform on [0, 1−2α], so the median 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(97)76068-3). diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index b592e9dc5..937391670 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -47,7 +47,7 @@ The methods draw on, and in places reimplement, solutions from: - P. Gasparotto, L. Barba, H.-C. Stadler et al., "TORO Indexer: a PyTorch-based indexing algorithm for kilohertz serial crystallography", *J. Appl. Cryst.* **57** (2024), 931-944 (the algorithm of the `ffbidx` fast-feedback indexer, §4). - I. Křivý & B. Gruber, "A unified algorithm for determining the reduced (Niggli) cell", *Acta Cryst.* **A32** (1976), 297-298, and International Tables for Crystallography Vol. A, Table 9.2.5.1 (the Niggli reduction and the lattice-character table of §5.3/§6). - J. E. Padilla & T. O. Yeates, "A statistic for local intensity differences: robustness to anisotropy and pseudo-centering and utility for detecting twinning", *Acta Cryst.* **D59** (2003), 1124-1130 (the L test, §13.2). -- T. O. Yeates, "Detecting and overcoming crystal twinning", *Methods Enzymol.* **276** (1997), 344-358 (the $H$ statistic behind the twin fraction of the twin hypothesis, §13.2). +- T. O. Yeates, "Detecting and overcoming crystal twinning", *Methods Enzymol.* **276** (1997), 344-358 (the twinned intensity model and the $H$ statistic behind the twin fraction of the twin hypothesis, §13.2, §14.9). - A. J. C. Wilson, "The probability distribution of X-ray intensities", *Acta Cryst.* **2** (1949), 318-321, and E. R. Howells, D. C. Phillips & D. Rogers, "The probability distribution of X-ray intensities. II. Experimental investigation and the X-ray detection of centres of symmetry", *Acta Cryst.* **3** (1950), 210-214 (the acentric and centric intensity distributions behind the centre-of-symmetry statistics of §13.1 and the Wilson outlier test of §13.3). - W. H. Baur & D. Kassner, "The perils of Cc: comparing the frequencies of falsely assigned space groups with their general population", *Acta Cryst.* **B48** (1992), 356-369 (writing the centrosymmetric group where the absences cannot decide the centre, §13.1). - R. J. Read, P. D. Adams & A. J. McCoy, "Intensity statistics in the presence of translational noncrystallographic symmetry", *Acta Cryst.* **D69** (2013), 176-183 (the native-Patterson detection of translational pseudo-symmetry, and the intensity modulation it produces, which the axial-zone screw-absence test scores against). diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index c6812bf84..ae3571ea6 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -64,7 +64,7 @@ angle < 0.44$) and the averaging closed a third or more of the gap from the unme A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500 and 0.333 untwinned, 0.375 and 0.200 for a perfect twin) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (2.0 for untwinned acentric data, 1.5 for a perfect twin) are written to the merged mmCIF as a twinning diagnostic. Both are taken on intensities divided by their resolution-shell mean, over the shells whose $\langle I/\sigma\rangle$ reaches 1, with Wilson outliers rejected; the $L$-test pairs are therefore compared on one scale even where two index steps span a steep fall-off, as in phenix.xtriage and ctruncate, and the selection is by shell, never by the individual reflection's $I/\sigma$, which would cut the weak tail and bias $\langle|L|\rangle$ down. The twin fraction is quoted from the statistic that carries the verdict — the $L$-test unless the call rests on the second moment alone — and a second moment is not turned into a fraction under a detected pseudo-translation, which inflates it. A twin law exists only where the **lattice** has a rotation the group lacks, so whether one is possible is read off the cell, not off the name of the Laue class: the lattice's two-folds are found as for the free set (Le Page, within 3° obliquity — GEMMI's `find_twin_laws`), and where none lies outside the group no twin is called. That covers the holohedral classes ($4/mmm$ on a tetragonal lattice, $6/mmm$, $m\bar{3}m$, $\bar{3}m$ on a rhombohedral one) and equally $mmm$ or $2/m$ on a lattice with no special metric, while a pseudo-merohedral $2/m$ or $mmm$ lattice — a monoclinic $\beta$ near 90°, a C-centred cell with $b \approx a\sqrt{3}$ — admits a law and stays eligible. The $L$-test is still read where no law is possible, for a different question: merging $I(h)$ with $I(Th)$ under a false operator $T$ gives $(I(h)+I(Th))/2$ whatever the twin fraction, which has the perfect-twin distribution, so $\langle|L|\rangle$ below 0.42 there is reported as **an adopted operator averaging unequal intensities** — the space group is too high, or a twin law was absorbed into the point group — a warning, never a change to the space group. A genuine operator leaves the untwinned 0.5. Reflections that overlap along a very long axis narrow the distribution the same way, which is one reason it stays a warning. The same numbers measured on the P1 merge the space-group search was given, before any point group was adopted, are reported beside them (`_BEFORE_SEARCH`), with that merge's own pseudo-translation declared and the reflections its lattice centring extinguishes left out. Between that bar and the twin band's 0.44 the verdict is `NO_TWIN_LAW`: the intensities read twin-like, but the lattice admits no law for the group, so they are not a twin in it — the group may be too high (a subgroup's twin law absorbed, or non-crystallographic symmetry taken for real), or something other than twinning narrowed the distribution. It is reported as such rather than as a twin that "needs a twin law" no refinement program could be given. -**The twin hypothesis.** Where the adopted group may have absorbed a twin law — a twin gate of the search refused it and was overruled (§13.1), or a group above $P1$ was adopted on a merge that read twin-like before the search — the same observations are also merged in a subgroup and written beside the answer as `_twin_.mtz` / `.cif`, in the subgroup's standard setting, to the answer's resolution and with the answer's free set (each reflection takes the flag of the reflection of the answer's merge its orbit falls in, so twin mates share it). The subgroup is the one a decision named — the group a refused promotion was adopted over, or the indexed lattice's group a metric re-ask promoted from — and otherwise, among the index-2 subgroups of the adopted group that the search confirmed as point-group finalists, the one whose lost two-fold agrees worst with the intensities (lowest operator CC): the most likely twin law. One file is written; the other subgroups that qualified are named in the report. The law is in the mmCIF's `_pdbx_reflns_twin` with its fraction, from Yeates' $H$ on the pairs it relates in the subgroup merge — $H = |I(h)-I(Th)|/(I(h)+I(Th))$ is uniform on $[0, 1-2\alpha]$, so $\alpha = 1/2 - \mathrm{median}\,H$, a lower bound since measurement error only adds to $H$. Nothing the answer is made of depends on it: it is a second merge from the same scaled observations, and the answer's files are byte for byte what they are without it. Refinement programs want twinned data merged in the true lower group, never detwinned, which is why the file is the merge and not a correction of the answer; a near-perfect twin of a subgroup and the higher group predict the same intensities, so which of the two is the crystal is a question for a model. +**The twin hypothesis.** Where the adopted group may have absorbed a twin law — a twin gate of the search refused it and was overruled (§13.1), or a group above $P1$ was adopted on a merge that read twin-like before the search — the same observations are also merged in a subgroup and written beside the answer as `_twin_.mtz` / `.cif`, in the subgroup's standard setting, to the answer's resolution and with the answer's free set (each reflection takes the flag of the reflection of the answer's merge its orbit falls in, so twin mates share it). The subgroup is the one a decision named — the group a refused promotion was adopted over, or the indexed lattice's group a metric re-ask promoted from — and otherwise, among the index-2 subgroups of the adopted group that the search confirmed as point-group finalists, the one whose lost two-fold agrees worst with the intensities (lowest operator CC): the most likely twin law. One file is written; the other subgroups that qualified are named in the report. The law is in the mmCIF's `_pdbx_reflns_twin` with its fraction, from Yeates' $H$ on the pairs it relates in the subgroup merge — $H = |I(h)-I(Th)|/(I(h)+I(Th))$ is uniform on $[0, 1-2\alpha]$, so $\alpha = 1/2 - \mathrm{median}\,H$, a lower bound since measurement error only adds to $H$. Nothing the answer is made of depends on it: it is a second merge from the same scaled observations, and the answer's files are byte for byte what they are without it. Refinement programs want twinned data merged in the true lower group, never detwinned, which is why the file is the merge and not a correction of the answer; a near-perfect twin of a subgroup and the higher group predict the same intensities, so which of the two is the crystal is a question for a model (§14.9). Beside them, on rotation data with `--finalist-ledger`, the run reports **twin-immune zone evidence** for the operators the adopted point group adds over each of its index-2 subgroups, read on a full-resolution merge of the observations in P1 that the option makes for the purpose. Under a twin law $T$, $I_\mathrm{obs}(h)=(1-\alpha)I(h)+\alpha I(Th)$; a reflection whose twin mate is itself up to the subgroup and Friedel is untouched at every $\alpha$, and those are exactly the reflections centric in the group but acentric in the subgroup. They read centric ($\langle|E^2-1|\rangle = 0.968$) if the added operators are real and acentric (0.736) if they are a twin law or a pseudo-symmetry — the one intensity statistic that still separates the two at $\alpha = 0.5$, where every operator statistic reads "real". Each zone is normalised against its own mean in resolution bins (and within the two phase classes of a detected pseudo-translation), read only in shells with $\langle I/\sigma\rangle \ge 5$ because noise inflates every class towards centric, and reported with $n$, a standard error and the centric-over-acentric Wilson log-likelihood ratio in nats (both densities convolved with each reflection's measurement error), beside the acentric control. It is read absolutely, never as the difference to the control: a perfect twin's control (0.541) plus noise inflation would otherwise look like true symmetry. An absolute reading is only as good as the normalisation, though, and the acentric control certifies it: an acentric population reads $-0.130$ nats per reflection when the normalisation is right, and whatever the control reads above that is the normalisation's — anisotropy, a pseudo-translation, a pseudo-centring, noise all inflate every class towards centric alike — so the zone, normalised the same way, carries the same per reflection and the **calibrated** evidence has it taken off (a twinned control reads *below* the expectation, so nothing is taken off a twin). The centric side is read **twinned at the fraction the lattice's other operations show** by their own correlation — the strongest CC of a lattice rotation outside the group, relative to the mean CC of the group's own operators, inverted through $\rho = 2\alpha(1-\alpha)/((1-\alpha)^2+\alpha^2)$ — because a twin by another law reaches a genuine zone exactly as it reaches the control: a genuine 321 crystal twinned by a 622 law read its 2-folds at −206 nats against the untwinned centric density and +209 at that law's fraction of 0.07. The operators' own twin law cannot reach their zone, so the acentric side stays untwinned, and the control's expectation is taken at the same fraction. Measured, a 6/m crystal with a 67 Ų anisotropy read its control at $+0.10$ nats per reflection and the zone of a refused 622 the same, so the zone's $+124$ nats were the normalisation's; calibrated it reads $-115$, and genuine promotions keep $+200$ and above. The calibrated zone is what a refused promotion is decided on (above). It is an in-house method: published practice (Yeates' $H$-test, xtriage) excludes these reflections as uninformative about the twin fraction, and no published test was found that uses them positively; the one precedent here is a trigonal crystal whose twin-immune zone read centric and whose higher group was then confirmed by refinement. @@ -273,3 +273,9 @@ The rotation and translation actually taken are reported (`RIGID_BODY_ROTATION_D The coordinates **as placed** are written as `_model.cif`, and as `_model.pdb` where the PDB format can hold the cell (the fragment-screening tools the file feeds — PanDDA, and the pair dimple produces — take a PDB beside the MTZ), so there is a coordinate file that describes the maps: the input's chains, residues, ligands, waters, B-factors, occupancies and anisotropic $U$s, moved, in the same cell and space group as `.mtz`. Taking the frame from the written reflections rather than from the input model matters — §14.6 may have relabelled them to the model's enantiomorph, which is neither the data's original group nor, necessarily, the model's. It is written **whenever the maps are**, not only where the rigid-body step was committed: the model is re-fractionalized into the data cell and may be relabelled whatever the placement decided, so an unmoved model is still not the input file; and a model the null of §14.5 rejected is scored, placed and mapped like any other, which is precisely the case where the density is worth looking at. +### 14.9 The twin hypothesis: which merge the maps come from + +Where the run wrote a twin hypothesis (§13.2), the model is scored against that merge too: as merged, and **as twinned** by its law — $|F_\mathrm{twin}(h)|^2 = (1-\alpha)|F_\mathrm{model}(h)|^2 + \alpha|F_\mathrm{model}(Th)|^2$ from the placed, scaled model, one overall scale refitted at each $\alpha$, $\alpha$ scanned from 0 to 0.5 in steps of 0.01 on the working set and R-free read at the chosen one. The placement, the anisotropic $B$ and the bulk solvent are those of the untwinned fit (§14.2, §14.8), so this is the same few-parameter comparison as §14.3, not a twin refinement. Against the answer's merge the model is scored as in §14.3, the subgroup model's structure factors simply read at each of the higher group's reflections. + +The maps and the placed model are then made from whichever merge the model fits better on R-free — the lower of the subgroup's plain and twinned R-free against the answer's — and the report says which (`MAPS_FROM=`) and with what R-free for each. The map coefficients are those of §14.4 against the chosen merge's amplitudes, not detwinned. **This chooses the maps and nothing else**: the space group, the answer's reflection files and every number in the report outside the twin keys are what the run decided without the model, because a model can be wrong and the reflection files must not depend on one (§14.5). The structure factors come from the same GPU engine as §14.1 where there is a GPU. + diff --git a/docs/RUGNUX_ADVANCED.md b/docs/RUGNUX_ADVANCED.md index 72a683ba6..0ac0a01f1 100644 --- a/docs/RUGNUX_ADVANCED.md +++ b/docs/RUGNUX_ADVANCED.md @@ -170,7 +170,15 @@ PDB: PanDDA's per-dataset input is `.pdb` beside `.mtz`, and dimple pair. (The PDB is skipped, with a log line, for a cell its fixed-width format cannot hold; the mmCIF is unconditional.) That is the file to open with the maps, and the one to hand to REFMAC5 or `phenix.refine`; its starting R-free is the `R_FREE=` the report quotes. -(`.cif` is the merged **reflections** — hence the separate name.) Both are written whenever the +(`.cif` is the merged **reflections** — hence the separate name.) + +Where the run also wrote a **twin hypothesis** — `_twin_.mtz`/`.cif`, the same data merged +in a subgroup with a twin law, written when the adopted space group may have absorbed one (see +[the twin hypothesis](RUGNUX_REPORT.md#the-twin-hypothesis-twin_subgroup_)) — the model is scored +against that merge as well, plainly and as twinned by the law, and the maps and `_model.cif` +come from whichever of the two merges the model fits better on R-free. The report says which +(`MAPS_FROM=`) and gives the R-free of both. The model only chooses the maps: `.mtz` and the +space group are the run's own either way. Both are written whenever the maps are, including for a model the data rejected: the rejection is a result, and it is exactly the case where someone wants to look at the model in the density. diff --git a/docs/RUGNUX_REPORT.md b/docs/RUGNUX_REPORT.md index 8634ca78b..fd76f8b88 100644 --- a/docs/RUGNUX_REPORT.md +++ b/docs/RUGNUX_REPORT.md @@ -861,6 +861,14 @@ as usual. | `MODEL_CHANGE_OF_BASIS`, `MODEL_SETTING_AS_READ` | Present only where the model was written in another description of the lattice (other axes, other centring, or another point group) and was put into the data's to be scored. Where such a model fits, the data are then written in the model's setting instead (`SETTING_SOURCE=MODEL`) and these keys are absent, since the model no longer moves. An alternative indexing of the same point group is settled here only where a reference MTZ (`-z`) has already fixed the data's indexing; otherwise the data are reindexed into the model's, and that is `MODEL_INDEXING_OPERATOR` | | `MODEL_INDEXING_MARGIN`, `MODEL_INDEXING_MARGIN_NULL`, `MODEL_INDEXING_MARGIN_SIGMA` | Present only where a merohedral ambiguity was probed. The winner's lead over the runner-up in R-free, against the lead a random placement of the same model produces | +Where a twin hypothesis was written (above), the model is scored against that merge too - as merged, +and as twinned by its law at a fraction scanned on the working set - and **the maps and the placed model +come from whichever merge the model fits better on R-free**. `MAPS_FROM=` says which (`MERGE` or +`TWIN_SUBGROUP_MERGE`), beside `TWIN_SUBGROUP_R_WORK`, `TWIN_SUBGROUP_R_FREE`, and, as twinned, +`TWIN_SUBGROUP_R_WORK_TWINNED`, `TWIN_SUBGROUP_R_FREE_TWINNED` and `TWIN_SUBGROUP_MODEL_TWIN_FRACTION`. +`R_WORK`/`R_FREE` above stay the answer's merge's. The model chooses the maps only: it never changes +the space group or `.mtz`. + 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/docs/RUGNUX_TUTORIAL.md b/docs/RUGNUX_TUTORIAL.md index 89eac3cfe..67988771c 100644 --- a/docs/RUGNUX_TUTORIAL.md +++ b/docs/RUGNUX_TUTORIAL.md @@ -287,7 +287,8 @@ Absorption is corrected only by the empirical surface a rotation run fits and ke [Validating against a model](RUGNUX_ADVANCED.md#validating-against-a-model-rugnux---model). - `_twin_.mtz` / `.cif` — only where the chosen space group may have absorbed a twin law: the same data merged in a subgroup ``, with the twin law and its fraction, same resolution and - free set as `.mtz`, for refinement with that law. See + free set as `.mtz`, for refinement with that law. With `--model`, the maps come from + whichever of the two merges the model fits better. See [The twin hypothesis](RUGNUX_REPORT.md#the-twin-hypothesis-twin_subgroup_). - `_report.txt` — the **results report**: what the run determined, in a form both a person and a beamline script can read. Always written, next to the files above. See diff --git a/image_analysis/structure_refinement/ModelValidation.cpp b/image_analysis/structure_refinement/ModelValidation.cpp index e355be145..ab8797727 100644 --- a/image_analysis/structure_refinement/ModelValidation.cpp +++ b/image_analysis/structure_refinement/ModelValidation.cpp @@ -423,6 +423,7 @@ ModelValidationResult Validate(const std::vector &merged, double wavelength_A, const std::vector &report_shell_d_min, const ModelValidationSchedule &schedule, + const gemmi::Op *twin_law, bool rigid_body_gpu) { ModelValidationResult result; result.model_path = model_path; @@ -1228,6 +1229,61 @@ ModelValidationResult Validate(const std::vector &merged, result.b_sol = best.b_sol; result.k_overall = best.k_overall; + // --- the twin target, where a twin law was given --- + // The intensity a twinned crystal records at h is (1-a) I(h) + a I(Th), so the model is compared + // against the data as twinned: |F_twin(h)|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 from the + // placed, scaled model, with one overall scale refitted at each a. The fraction is scanned in steps + // of 0.01 and chosen on the working set; R-free is read at it. Nothing else is refitted against the + // twinned target - the placement, the anisotropic B and the bulk solvent are those of the untwinned + // fit - so this is the same few-parameter comparison as R above, not a twin refinement. + // Following Yeates (1997) Methods Enzymol. 276, 344-358 + if (twin_law != nullptr) { + std::unordered_map model_i; + model_i.reserve(fmodel.v.size()); + for (const auto &hv : fmodel.v) + model_i[hkl_key(hv.hkl)] = std::norm(hv.value); + struct TwinTerm { double fo, i1, i2; bool free; }; + std::vector twin_terms; + for (const auto &hv : fmodel.v) { + const auto it = obs_by_hkl.find(hkl_key(hv.hkl)); + if (it == obs_by_hkl.end()) continue; + const gemmi::Miller mate = asu.to_asu(twin_law->apply_to_hkl(hv.hkl), gops).first; + const auto m = model_i.find(hkl_key(mate)); + if (m == model_i.end()) continue; + twin_terms.push_back({it->second.F, std::norm(hv.value), m->second, it->second.free}); + } + double best_r_work = 2.0; + for (int step = 0; step <= 50 && !twin_terms.empty(); ++step) { + const double a = 0.01 * step; + double fo_fc = 0, fc_fc = 0; + for (const auto &t : twin_terms) + if (!t.free) { + const double fc = std::sqrt((1.0 - a) * t.i1 + a * t.i2); + fo_fc += t.fo * fc; + fc_fc += fc * fc; + } + const double k = fc_fc > 0 ? fo_fc / fc_fc : 1.0; + double num_w = 0, den_w = 0, num_f = 0, den_f = 0; + for (const auto &t : twin_terms) { + const double diff = std::fabs(t.fo - k * std::sqrt((1.0 - a) * t.i1 + a * t.i2)); + if (t.free) { num_f += diff; den_f += t.fo; } + else { num_w += diff; den_w += t.fo; } + } + if (den_w > 0 && den_f > 0 && num_w / den_w < best_r_work) { + best_r_work = num_w / den_w; + result.r_work_twin = best_r_work; + result.r_free_twin = num_f / den_f; + result.twin_fraction = a; + } + } + result.twin_target = std::isfinite(result.r_free_twin); + if (result.twin_target) + logger.Info("Model validation: against the data as twinned by {} - fraction {:.2f} (scanned on the " + "working set) - R-work {:.4f}, R-free {:.4f}; untwinned {:.4f} / {:.4f}", + twin_law->as_hkl().triplet('h'), result.twin_fraction, result.r_work_twin, + result.r_free_twin, result.r_work, result.r_free); + } + // --- CC(model, data) by resolution shell --- // The correlation of the merged intensities with |F_model|^2, binned on the merge table's own // shells so the two tables line up row for row. Nearly free: the scaled Fmodel is already here, @@ -1655,14 +1711,16 @@ ModelValidationResult ValidateAgainstModel(const std::vector & size_t nthreads, double wavelength_A, const std::vector &report_shell_d_min, - const ModelValidationSchedule &schedule) { + const ModelValidationSchedule &schedule, + const gemmi::Op *twin_law) { #ifdef JFJOCH_USE_CUDA // A CUDA failure ends the validation, not the run, and nothing is moved to the CPU part way through: // the validation reports why it did not finish, and a validation that did not finish decides nothing, // so the reflection files are those of a run without a model. try { return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, - probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, true); + probe_indexing_ambiguity, nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, + true); } catch (const RigidBodyGPUFailure &e) { cuda_clear_error(); ModelValidationResult failed; @@ -1673,7 +1731,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & } #else return Validate(merged, cell, model_path, output_prefix, logger, data_space_group, probe_indexing_ambiguity, - nthreads, wavelength_A, report_shell_d_min, schedule, false); + nthreads, wavelength_A, report_shell_d_min, schedule, twin_law, false); #endif } diff --git a/image_analysis/structure_refinement/ModelValidation.h b/image_analysis/structure_refinement/ModelValidation.h index c89feb567..adc686bf2 100644 --- a/image_analysis/structure_refinement/ModelValidation.h +++ b/image_analysis/structure_refinement/ModelValidation.h @@ -184,6 +184,14 @@ struct ModelValidationResult { bool adopted_model_enantiomorph = false; int model_space_group_number = 0; + // Against the data as twinned by the twin law ValidateAgainstModel was given (none given: not + // made): R-work and R-free with |F_twin|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2 and one overall + // scale, at the fraction a (0 to 0.5) that gives the lowest R-work. + bool twin_target = false; + double twin_fraction = NAN; + double r_work_twin = NAN; + double r_free_twin = NAN; + // The model as it was scored: re-fractionalized into the data cell and moved by the rigid-body // step, everything else exactly as it was read. Null where the validation did not run. The frame // it is written in is only settled by AdoptModelFrame below, so it is the caller that puts it on @@ -224,6 +232,9 @@ struct ModelValidationResult { // nthreads is what the null's replicates and the candidate changes of basis run on - each is // independent of the others, so they run at once. The answer does not depend on it. // +// twin_law, where given, is a twin law of the data's group (a Miller-index operator in the data's +// basis): the model is then also scored against the data as twinned by it (twin_target and below). +// // report_shell_d_min is the merge statistics' own shell bounds, coarse to fine, and it is what // cc_model_shells is binned on. Sharing the grid is the point: a reader has to be able to put a // CC(model, data) row beside that shell's CC1/2 and know the two describe the same reflections. @@ -261,7 +272,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & size_t nthreads = 1, double wavelength_A = 0.0, const std::vector &report_shell_d_min = {}, - const ModelValidationSchedule &schedule = {}); + const ModelValidationSchedule &schedule = {}, + const gemmi::Op *twin_law = nullptr); // Reindex `merged` into the frame ValidateAgainstModel reported, so the reflection files that are // written describe the same indexing as the R-factors and the maps. Returns the space group they are diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 497d3150b..d8f50ed55 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -2392,6 +2392,33 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, } if (!mv.maps_prefix.empty()) Add(s, KeyText("MAPS_PREFIX", mv.maps_prefix)); + // 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. + if (result.twin_subgroup_merge && result.twin_subgroup_merge->model_validation) { + const auto &t = *result.twin_subgroup_merge; + const auto &tv = *t.model_validation; + Add(s, KeyEnum("MAPS_FROM", t.maps_from_this_merge ? "TWIN_SUBGROUP_MERGE" : "MERGE", + {"MERGE", "TWIN_SUBGROUP_MERGE"})); + if (tv.ok) { + Add(s, KeyReal("TWIN_SUBGROUP_R_WORK", tv.r_work, "{:.4f}")); + Add(s, KeyReal("TWIN_SUBGROUP_R_FREE", tv.r_free, "{:.4f}")); + if (tv.twin_target) { + Add(s, KeyReal("TWIN_SUBGROUP_R_WORK_TWINNED", tv.r_work_twin, "{:.4f}")); + Add(s, KeyReal("TWIN_SUBGROUP_R_FREE_TWINNED", tv.r_free_twin, "{:.4f}")); + Add(s, KeyReal("TWIN_SUBGROUP_MODEL_TWIN_FRACTION", tv.twin_fraction, "{:.2f}")); + } + } + Add(s, Blank()); + Add(s, Prose(fmt::format( + " The model was also scored against the twin hypothesis' merge in {} ({}.mtz): as merged,\n" + " and as twinned by {} - |F_twin|^2 = (1-a)|F_model(h)|^2 + a|F_model(Th)|^2, the\n" + " fraction a scanned on the working set. The maps and the placed model come from {},\n" + " the merge the model fits better on R-free ({:.4f} against {:.4f}). This chooses the\n" + " maps only: the answer's reflection files and every number above are the answer's.", + t.space_group, t.file_prefix, t.twin_law, + t.maps_from_this_merge ? "the merge in " + t.space_group : "the answer's merge", + tv.twin_target ? std::min(tv.r_free, tv.r_free_twin) : tv.r_free, mv.r_free))); + } if (!mv.cc_model_shells.empty()) { Add(s, Blank()); diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index c6eb42c66..782347db0 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -491,6 +491,10 @@ struct ProcessResult { double twin_fraction = NAN; // EstimateTwinLawFraction: a lower bound int twin_fraction_pairs = 0; std::vector other_subgroups; // other confirmed subgroups that qualified + // --model: the model scored against this merge (plain and as twinned by the law), and whether + // the maps and the placed model come from it rather than from the answer's merge. + std::optional model_validation; + bool maps_from_this_merge = false; }; std::optional twin_subgroup_merge; diff --git a/rugnux/RugnuxScaleMerge.cpp b/rugnux/RugnuxScaleMerge.cpp index 12641df51..4dde9ed34 100644 --- a/rugnux/RugnuxScaleMerge.cpp +++ b/rugnux/RugnuxScaleMerge.cpp @@ -3738,7 +3738,61 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) { // file". A rejected model is scored, placed and mapped like any other, and that is the case // where someone most wants to see the model in the density, so it gets its file too. // Written beside the reflection files below: nothing after this reads it. - if (validation.placed_model) + // The model asked about the twin hypothesis too, where one was written: scored against the + // subgroup's merge, plainly and as twinned by its law (the twin target, at a scanned fraction), + // against the answer's merge as above. The maps and the placed model come from whichever merge + // the model fits better on R-free; the answer's files are the answer's either way - the model + // chooses the maps, never the merge (see docs/RUGNUX_REPORT.md). + prepare_twin_file(); + bool maps_from_twin = false; + if (twin_file && validation.ok) { + const double answer_r_free = validation.r_free; + const auto twin_r_free = [](const ModelValidationResult &v) { + return v.twin_target ? std::min(v.r_free, v.r_free_twin) : v.r_free; + }; + ModelValidationSchedule twin_schedule; + twin_schedule.maps_superseded = [&](const ModelValidationResult &v) { + return !(twin_r_free(v) < answer_r_free); + }; + phase("Validating against model (twin hypothesis)"); + auto twin_validation = ValidateAgainstModel( + twin_file->merged, twin_file->cell, config_.model_path, config_.output_prefix, logger, + &twin_file->space_group, /*probe_indexing_ambiguity=*/config_.reference_data.empty(), + static_cast(config_.nthreads), experiment_.GetWavelength_A(), {}, twin_schedule, + &twin_file->law); + maps_from_twin = twin_validation.ok && twin_r_free(twin_validation) < answer_r_free; + // The subgroup file follows its own validation's indexing, as the answer follows its own: + // the file, the R-factors and the maps then describe one indexing. + if (twin_validation.ok) { + const gemmi::Op op = twin_validation.indexing_op; + twin_file->space_group = *AdoptModelFrame(twin_validation, twin_file->merged, + twin_file->space_group, + experiment_.GetScalingSettings().GetMergeFriedel(), + logger); + if (!(op == gemmi::Op::identity())) + twin_file->law = op.inverse().combine(twin_file->law).combine(op); + result.twin_subgroup_merge->twin_law = twin_file->law.as_hkl().triplet(); + result.twin_subgroup_merge->space_group = twin_file->space_group.xhm(); + } + logger.Info("Model validation, twin hypothesis: R-free {:.4f} against the answer's merge in {}, " + "{:.4f} against the merge in {}{} - the maps and the placed model come from {}", + answer_r_free, experiment_.GetSpaceGroupOrP1().xhm(), twin_validation.r_free, + twin_file->space_group.xhm(), + twin_validation.twin_target + ? fmt::format(" and {:.4f} as twinned at a fraction of {:.2f}", + twin_validation.r_free_twin, twin_validation.twin_fraction) + : std::string(), + maps_from_twin ? "the merge in " + twin_file->space_group.xhm() : "the answer's merge"); + if (maps_from_twin && twin_validation.placed_model) + placed_model_file = std::async(std::launch::async, + [placed = twin_validation.placed_model, cell = twin_file->cell, + sg = twin_file->space_group, prefix = config_.output_prefix, &logger] { + WritePlacedModel(*placed, cell, sg, prefix, logger); + }); + result.twin_subgroup_merge->maps_from_this_merge = maps_from_twin; + result.twin_subgroup_merge->model_validation = std::move(twin_validation); + } + if (validation.placed_model && !maps_from_twin) placed_model_file = std::async(std::launch::async, [placed = validation.placed_model, cell = *result.consensus_cell, sg = experiment_.GetSpaceGroupOrP1(), prefix = config_.output_prefix, &logger] {