From a4cc908ff1ff38acf64790ae4a675f7c31ec367b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 9 Jul 2026 21:59:37 +0200 Subject: [PATCH] rugnux: run the L-test on rhombohedral data + write twinning stats to CIF The Padilla-Yeates L-test paired each reflection with a neighbour two steps along an axis ((h+2,k,l) etc.). A step of 2 preserves the reflection condition of P/I/C/F/A/B lattices but violates R-centring (-h+k+l = 0 mod 3, since 2 != 0 mod 3), so every neighbour was systematically absent -> zero pairs, and the L-test silently did not run on ANY rhombohedral crystal (twinning there was left to the weaker second-moment test alone). Add diagonal (1,1,0)/(1,1,3) fallback steps that preserve the mod-3 condition in both obverse and reverse settings, tried only when the axis steps find no present partner (P/I/C/F behaviour unchanged). On the twinned R3 insulin the L-test now runs: 25k pairs, = 0.224 -> twinning detected. Also write the twinning indicators to the output mmCIF next to ISa: _reflns.pdbx_L_test_mean_abs_L / _L_test_mean_L_squared / _second_moment_I, as free-text pdbx values, guarded on pair/reflection count. Threaded the TwinningAnalysisResult into WriteReflections (both the full-pipeline and --scale callers). Note: on rhombohedral lattices the L-test reads ~0.03 lower than on P/I/C/F (its allowed step is slightly more local), so it separates the insulin twin pair only weakly (twinned 0.398 vs untwinned 0.427) and the second moment remains the cleaner discriminator there. Co-Authored-By: Claude Opus 4.8 (1M context) --- image_analysis/WriteReflections.cpp | 15 ++++++++++++++- image_analysis/WriteReflections.h | 3 +++ image_analysis/scale_merge/TwinningAnalysis.cpp | 9 ++++++++- rugnux/Rugnux.cpp | 2 +- tools/rugnux_cli.cpp | 5 +++-- 5 files changed, 29 insertions(+), 5 deletions(-) diff --git a/image_analysis/WriteReflections.cpp b/image_analysis/WriteReflections.cpp index 826d5ba7..15348a57 100644 --- a/image_analysis/WriteReflections.cpp +++ b/image_analysis/WriteReflections.cpp @@ -3,6 +3,7 @@ #include "WriteReflections.h" #include "scale_merge/Merge.h" +#include "scale_merge/TwinningAnalysis.h" #include #include @@ -62,6 +63,7 @@ void WriteMmcifReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename) { std::ofstream out(filename); @@ -136,6 +138,16 @@ void WriteMmcifReflections(const std::vector &reflections, out << "_reflns.pdbx_Rrim_I_all " << Fmt(ov.r_meas, 4) << "\n"; out << "_reflns.pdbx_CC_half " << Fmt(ov.cc_half, 4) << "\n"; out << "_reflns.pdbx_diffrn_ISa " << CifStr(isa) << " # asymptotic I/sigma (Diederichs)\n"; + // Twinning indicators (no standard CIF item, written as free-text pdbx values like ISa above). + if (twinning.l_test_pairs > 0) { + out << "_reflns.pdbx_L_test_mean_abs_L " << Fmt(twinning.mean_abs_l, 3) + << " # Padilla-Yeates <|L|> (untwinned 0.500, perfect twin 0.375)\n"; + out << "_reflns.pdbx_L_test_mean_L_squared " << Fmt(twinning.mean_l_squared, 3) + << " # (untwinned 0.333, perfect twin 0.200)\n"; + } + if (twinning.moment_reflections > 0) + out << "_reflns.pdbx_second_moment_I " << Fmt(twinning.second_moment, 3) + << " # /^2 (untwinned 2.00, perfect twin 1.50)\n"; out << "#\n"; out << "loop_\n"; @@ -249,13 +261,14 @@ void WriteReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename) { switch (experiment.GetScalingSettings().GetFileFormat()) { case IntensityFormat::Text: WriteHKLReflections(reflections, filename + ".hkl"); break; case IntensityFormat::mmCIF: - WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, filename + ".cif"); + WriteMmcifReflections(reflections, unitCell, experiment, statistics, isa, twinning, filename + ".cif"); break; case IntensityFormat::MTZ: WriteMtzReflections(reflections, unitCell, experiment, filename + ".mtz"); diff --git a/image_analysis/WriteReflections.h b/image_analysis/WriteReflections.h index 8215dd93..4c03d092 100644 --- a/image_analysis/WriteReflections.h +++ b/image_analysis/WriteReflections.h @@ -11,12 +11,14 @@ #include "../common/DiffractionExperiment.h" struct MergeStatistics; +struct TwinningAnalysisResult; void WriteMmcifReflections(const std::vector &reflections, const UnitCell &unitCell, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename); void WriteMtzReflections(const std::vector &reflections, @@ -32,4 +34,5 @@ void WriteReflections(const std::vector &reflections, const DiffractionExperiment &experiment, const MergeStatistics &statistics, const std::string &isa, + const TwinningAnalysisResult &twinning, const std::string &filename); \ No newline at end of file diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 97b70c30..8dd392d7 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -49,7 +49,14 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector& merg if (UsableIntensity(r) && r.I > 0.0 && acentric(r)) intensity.emplace(PackHKL(r.h, r.k, r.l), r.I); - const std::array, 3> steps{{{2, 0, 0}, {0, 2, 0}, {0, 0, 2}}}; + // Step to a nearby, symmetry-independent partner. The axis step of 2 preserves the reflection + // condition of P/I/C/F/A/B lattices (parity-based), but violates R-centring (-h+k+l = 0 mod 3, + // as 2 != 0 mod 3) -> the partner is systematically absent and rhombohedral crystals yield zero + // pairs. The diagonal (1,1,0)/(1,1,3) steps preserve the mod-3 condition in BOTH obverse and + // reverse settings; they are tried only when the axis steps find no present partner, so P/I/C/F + // behaviour is unchanged (first present partner wins). + const std::array, 5> steps{ + {{2, 0, 0}, {0, 2, 0}, {0, 0, 2}, {1, 1, 0}, {1, 1, 3}}}; double sum_abs_l = 0.0, sum_l2 = 0.0; int n_pairs = 0; diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 80b5dd9c..fae65aa4 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -609,7 +609,7 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { phase("Writing reflections"); WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, result.error_model_isa > 0 ? fmt::format("{:.2f}", result.error_model_isa) : "?", - config_.output_prefix); + result.twinning, config_.output_prefix); // Per-image scaling table (G, B-factor, mosaicity, wedge, CC) for inspection / XDS // comparison. The offline self-scaling result is otherwise not exposed (process.h5's // per-image arrays are only filled on the online per-image path). Sourced from the diff --git a/tools/rugnux_cli.cpp b/tools/rugnux_cli.cpp index bb3563ad..49c12685 100644 --- a/tools/rugnux_cli.cpp +++ b/tools/rugnux_cli.cpp @@ -1014,12 +1014,13 @@ int main(int argc, char **argv) { const auto twin_sg_number = experiment.GetSpaceGroupNumber(); const gemmi::SpaceGroup *twin_sg = twin_sg_number ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; - std::cout << std::endl << TwinningAnalysisToText(AnalyzeTwinning(merged_reflections, twin_sg)) << std::endl; + const auto twinning = AnalyzeTwinning(merged_reflections, twin_sg); + std::cout << std::endl << TwinningAnalysisToText(twinning) << std::endl; if (!output_prefix.empty()) WriteReflections(merged_reflections, *experiment.GetUnitCell(), experiment, merged_statistics, error_model_isa > 0 ? fmt::format("{:.2f}", error_model_isa) : "?", - output_prefix); + twinning, output_prefix); return 0; }