From 21dd99524f2b3c5fde0357af55892b6824bbcc6a Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Tue, 30 Jun 2026 14:18:49 +0200
Subject: [PATCH] rotation_indexer: map FFT primitive cell to the SG
 conventional setting

When a space group is supplied without a reference cell, de-novo two-pass
rotation indexing fed the FFT's Niggli-reduced primitive cell straight into
XtalOptimizer as if it were the conventional cell. For non-primitive lattices
(centered I/F/R/C, or hexagonal where the primitive pair sits at gamma=60) the
conventional-system model then refined to a wrong minimum and indexed 0% of
frames: cytC (P3121) gave 103.9/103.9/78 instead of 83.7/83.7/88.6, insulin
(I213) 66.7 instead of 77.65, insulin-R3 51/51/36 instead of 81.4/81.4/33.3.

Run LatticeSearch on the FFT primitive cell (it already yields the correct
conventional cell + reindex for I/F/R/C). For the one remaining gap - a
metrically hexagonal lattice that the geometry-keyed search lands on the
ortho-hexagonal C setting - re-express the reduced primitive cell in
conventional hexagonal axes (b -> b - a opens gamma 60 -> 120).

De-novo "-S" now indexes cytC/insu/Ins_H/lyso/MyoB/EP/lyso_ref at 100% with the
correct cell; the "-C -S" path is unchanged. The helper stays in this .cpp
(g++) rather than CrystalLattice.h to avoid recompiling CUDA units, which is
broken under the box's CUDA-13 nvcc; promote it to a method once that is fixed.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 .../rotation_indexer/RotationIndexer.cpp      | 35 +++++++++++++++++--
 1 file changed, 32 insertions(+), 3 deletions(-)

diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp
index f7a6de3b..2b40f7e7 100644
--- a/image_analysis/rotation_indexer/RotationIndexer.cpp
+++ b/image_analysis/rotation_indexer/RotationIndexer.cpp
@@ -9,6 +9,19 @@
 #include "../lattice_search/LatticeSearch.h"
 #include "../indexing/MultiLatticeSearch.h"
 
+namespace {
+    // Re-express a primitive hexagonal/trigonal lattice in the conventional hexagonal setting
+    // (a = b, gamma = 120). The Niggli-reduced primitive cell carries the two equal-length axes
+    // at gamma = 60; replacing b with b - a opens that angle to 120 without changing the lattice.
+    CrystalLattice HexagonalConventional(CrystalLattice latt) {
+        latt.ReorderABEqual();                       // put the equal-length pair in a, b
+        Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2();
+        if (angle_deg(a, b) < 90.0f)
+            b -= a;
+        return CrystalLattice(a, b, c);              // constructor fixes handedness
+    }
+}
+
 RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer)
     : experiment(x),
       index_ice_rings(x.GetIndexingSettings().GetIndexIceRings()),
@@ -40,12 +53,28 @@ void RotationIndexer::RunIndexing() {
     if (!indexer_result.lattice.empty() && indexer_result.lattice[0].CalcVolume() > 1.0) {
         auto sg = experiment.GetGemmiSpaceGroup();
         if (sg) {
+            // The FFT returns a primitive cell; map it to the SG's conventional setting via the
+            // Niggli-class table (which handles I/F/R/C centering). LatticeSearch keys off geometry
+            // alone, so a metrically-hexagonal lattice can land on the ortho-hexagonal C setting -
+            // when the SG is trigonal/hexagonal, re-express the reduced primitive cell in the
+            // conventional hexagonal axes instead.
+            auto ls = LatticeSearch(indexer_result.lattice[0]);
+
+            const auto is_hexagonal = [](gemmi::CrystalSystem s) {
+                return s == gemmi::CrystalSystem::Trigonal || s == gemmi::CrystalSystem::Hexagonal;
+            };
+
+            CrystalLattice conventional = ls.conventional;
+            if (is_hexagonal(sg->crystal_system()) && !is_hexagonal(ls.system))
+                conventional = HexagonalConventional(ls.primitive_reduced);
+
             search_result_ = LatticeSearchResult{
-                .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one
-                .conventional = indexer_result.lattice[0], // If lattice provided, it is for now primitive == conventional
+                .niggli_class = ls.niggli_class,
+                .primitive_reduced = ls.primitive_reduced,
+                .conventional = conventional,
                 .system = sg->crystal_system(),
                 .centering = sg->centring_type(),
-                .reindex = gemmi::Mat33(1, 0, 0, 0, 1, 0, 0, 0, 1),
+                .reindex = ls.reindex,
             };
         } else {
             // Find lattice type based on cell