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