From fea2790a41163eb244493f4fc8ea407b1866f603 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 30 Jun 2026 21:31:20 +0200 Subject: [PATCH] rotation_indexer: detect hexagonal metric de-novo (fixes trigonal cytC indexing) The hexagonal lattice metric (two equal axes at 60/120 deg, both perpendicular to the third) is also satisfied by its ortho-hexagonal C-centred supercell, so the geometry-keyed LatticeSearch lands a metrically-hexagonal lattice on the C-orthorhombic setting. The -S path already re-expresses it (HexagonalConventional) keyed on the supplied space group; do the same de-novo, keyed on the metric of the reduced PRIMITIVE cell (rhombohedral lattices have a rhombohedral primitive cell and are unaffected). cytC (P3121) now indexes de-novo as hexagonal 83.8/83.8/88.6 gamma=120 and merges in P3121 (CC1/2 99.7%, ISa 13.1), for both test crystals, instead of the orthohexagonal C2/P1. lysoC/InsI3/InsH3 unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../rotation_indexer/RotationIndexer.cpp | 33 +++++++++++++++++-- 1 file changed, 31 insertions(+), 2 deletions(-) diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 2b40f7e7..4334abb9 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -20,6 +20,27 @@ namespace { b -= a; return CrystalLattice(a, b, c); // constructor fixes handedness } + + bool IsHexagonalSystem(gemmi::CrystalSystem s) { + return s == gemmi::CrystalSystem::Trigonal || s == gemmi::CrystalSystem::Hexagonal; + } + + // The hexagonal lattice metric (two equal axes at 60/120 deg, both perpendicular to the third) is + // also satisfied by its ortho-hexagonal C-centred supercell, so the geometry-keyed LatticeSearch can + // land there. Detect the hexagonal metric on the reduced PRIMITIVE cell so the de-novo path (no space + // group to key on) can re-express it in conventional hexagonal axes. + bool IsMetricallyHexagonal(CrystalLattice latt, float rel_tol = 0.03f, float angle_tol_deg = 3.0f) { + latt.ReorderABEqual(); + const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2(); + const float la = a.Length(), lb = b.Length(); + if (la <= 0.0f || lb <= 0.0f || std::fabs(la - lb) > rel_tol * std::max(la, lb)) + return false; + const float gab = angle_deg(a, b); + if (std::fabs(gab - 60.0f) > angle_tol_deg && std::fabs(gab - 120.0f) > angle_tol_deg) + return false; + return std::fabs(angle_deg(a, c) - 90.0f) <= angle_tol_deg && + std::fabs(angle_deg(b, c) - 90.0f) <= angle_tol_deg; + } } RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer) @@ -77,8 +98,16 @@ void RotationIndexer::RunIndexing() { .reindex = ls.reindex, }; } else { - // Find lattice type based on cell - search_result_ = LatticeSearch(indexer_result.lattice[0]); + // Find lattice type from the metric. LatticeSearch keys off geometry, so a metrically + // hexagonal lattice can land on the ortho-hexagonal C setting; re-express it in conventional + // hexagonal axes (P, gamma=120) so the trigonal/hexagonal 3-fold is not hidden from scaling. + auto ls = LatticeSearch(indexer_result.lattice[0]); + if (!IsHexagonalSystem(ls.system) && IsMetricallyHexagonal(ls.primitive_reduced)) { + ls.conventional = HexagonalConventional(ls.primitive_reduced); + ls.system = gemmi::CrystalSystem::Hexagonal; + ls.centering = 'P'; + } + search_result_ = ls; } // Run refinement