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, <L^2> = 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) <noreply@anthropic.com>

image_analysis/WriteReflections.cpp

@@ -3,6 +3,7 @@
 
 #include "WriteReflections.h"
 #include "scale_merge/Merge.h"
+#include "scale_merge/TwinningAnalysis.h"
 
 #include <cmath>
 #include <fstream>
@@ -62,6 +63,7 @@ void WriteMmcifReflections(const std::vector<MergedReflection> &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<MergedReflection> &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)
+                << "   # <L^2> (untwinned 0.333, perfect twin 0.200)\n";
+        }
+        if (twinning.moment_reflections > 0)
+            out << "_reflns.pdbx_second_moment_I       " << Fmt(twinning.second_moment, 3)
+                << "   # <I^2>/<I>^2 (untwinned 2.00, perfect twin 1.50)\n";
         out << "#\n";
 
         out << "loop_\n";
@@ -249,13 +261,14 @@ void WriteReflections(const std::vector<MergedReflection> &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");

image_analysis/WriteReflections.h

@@ -11,12 +11,14 @@
 #include "../common/DiffractionExperiment.h"
 
 struct MergeStatistics;
+struct TwinningAnalysisResult;
 
 void WriteMmcifReflections(const std::vector<MergedReflection> &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<MergedReflection> &reflections,
@@ -32,4 +34,5 @@ void WriteReflections(const std::vector<MergedReflection> &reflections,
                          const DiffractionExperiment &experiment,
                          const MergeStatistics &statistics,
                          const std::string &isa,
+                         const TwinningAnalysisResult &twinning,
                          const std::string &filename);

image_analysis/scale_merge/TwinningAnalysis.cpp

@@ -49,7 +49,14 @@ TwinningAnalysisResult AnalyzeTwinning(const std::vector<MergedReflection>& merg
         if (UsableIntensity(r) && r.I > 0.0 && acentric(r))
             intensity.emplace(PackHKL(r.h, r.k, r.l), r.I);
 
-    const std::array<std::array<int, 3>, 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<std::array<int, 3>, 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;

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

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;
     }