From 3c001eb7a3f50d6ac29f18070dbd4a489e4a187a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 11 Nov 2026 11:56:17 +0100 Subject: [PATCH] RESULT-CHANGING: LatticeSearch takes the Bravais class from the metric two-folds where the table misses it The 44-character table is read on the single Niggli-reduced cell. A lattice whose shortest vectors are degenerate has several reduced cells: an F-centred cubic one reduces to the 60/60/60 rhombohedral cell or to a 60/60/90 one (also the reduced cell of an I-tetragonal lattice with c = a*sqrt(2)) according to the sign of its strain, and only the first carries character 1. On an F-cubic small-molecule crystal the indexing therefore named the lattice I-tetragonal, the space-group search never got the cubic operators, -C with the cubic cell ended in P1 and -S Fm-3m was refused as "no setting of that lattice carries the group's". LatticeSearch now runs LePageLattice (basis-independent, already used by the space-group stage) on the same lattice and takes its class where it is richer than the table's and its conventional cell meets that class's metric within the caller's own length and angle tolerances (obliquity at most min(angle tolerance, 1 deg)). A tetragonal or cubic class refused only on its length equality is returned as the orthorhombic lattice of the same centring, so the demoted hypothesis (dist_tolerance 0) gets oF/oI/oP rather than triclinic. Measured on synthetic lattices, 300 per Bravais type, strain 0.05-0.5%: cF is now right 296-300/300 (was 155-159); hP 291-300 (105-108); oC, oI, tI, hR also gain, the monoclinic under-calls of centred orthorhombic lattices mostly go. Costs: a few more promotions within the existing tolerances (mP -> oC 39 -> 50 of 300), which the first pass's constrained-vs-free guard and the metric-symmetry check judge on the data. ci tier vs the rc184 full run: 29/29 pass, the null passes, 27 sets identical; myob_x10sa changes slightly (its first-pass candidate is pseudo-oC at 0.3 deg obliquity, now promoted and then dropped by the guard; P21, R_free 0.246 -> 0.245), 6iu5 d_min 2.12 -> 2.11. LatticePointGroupOrder moves from RotationIndexer to LatticeSearch; the class metric departure is shared with MetricDeviation. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- .../lattice_search/LatticeSearch.cpp | 99 +++++++++++++++---- image_analysis/lattice_search/LatticeSearch.h | 4 + .../rotation_indexer/RotationIndexer.cpp | 14 --- tests/LatticeSearchTest.cpp | 25 +++++ 4 files changed, 107 insertions(+), 35 deletions(-) diff --git a/image_analysis/lattice_search/LatticeSearch.cpp b/image_analysis/lattice_search/LatticeSearch.cpp index 7e5624a3f..b60c3c388 100644 --- a/image_analysis/lattice_search/LatticeSearch.cpp +++ b/image_analysis/lattice_search/LatticeSearch.cpp @@ -3,6 +3,7 @@ #include "../../common/JFJochMath.h" #include "LatticeSearch.h" +#include "LePageLattice.h" #include #include @@ -35,6 +36,27 @@ struct NiggliClass { }; namespace { +// How far a cell is from the metric its crystal system imposes: the largest angle departure (degrees) +// and the largest relative length difference among the axes the system makes equal. +std::pair ClassMetricDeparture(const UnitCell &uc, gemmi::CrystalSystem system) { + double angle = 0.0, length = 0.0; + auto a_dev = [&](double x, double ideal) { angle = std::max(angle, std::fabs(x - ideal)); }; + auto l_dev = [&](double x, double y) { length = std::max(length, std::fabs(x - y) / (0.5 * (x + y))); }; + switch (system) { + case gemmi::CrystalSystem::Monoclinic: a_dev(uc.alpha, 90); a_dev(uc.gamma, 90); break; + case gemmi::CrystalSystem::Orthorhombic: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 90); break; + case gemmi::CrystalSystem::Tetragonal: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 90); + l_dev(uc.a, uc.b); break; + case gemmi::CrystalSystem::Trigonal: + case gemmi::CrystalSystem::Hexagonal: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 120); + l_dev(uc.a, uc.b); break; + case gemmi::CrystalSystem::Cubic: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 90); + l_dev(uc.a, uc.b); l_dev(uc.b, uc.c); break; + default: break; // triclinic: nothing is asserted, nothing is violated + } + return {angle, length}; +} + // The body of both entry points. only_class, when given, keeps just the characters of that Bravais // class - see LatticeSearchForClass. With no filter this is the original walk unchanged, and the // triclinic character fits every metric, so it always returns a result. @@ -417,8 +439,6 @@ std::optional SearchCharacters(const CrystalLattice &L, dou // Character 44 fits any cell, so a match is always found - "nothing fits" is reported as triclinic. auto found = match(L_niggli, D, E, F); - if (found && found->system != gemmi::CrystalSystem::Triclinic) - return *found; // A reduced cell with an angle 90 to within NIGGLI_TYPE_BOUNDARY_DEG sits ON the boundary between the // two Niggli types: the same lattice reduces to an all-acute cell or an all-obtuse one according to @@ -433,7 +453,7 @@ std::optional SearchCharacters(const CrystalLattice &L, dou // run) and on an I-centred orthorhombic one whose reduced gamma sits 0.1 deg from 90. Beta is tried // first, so a cell the earlier beta-only retry already rescued is answered exactly as before. const UnitCell uc_niggli = L_niggli.GetUnitCell(); - if (D > 0 && E > 0 && F > 0) { + if ((!found || found->system == gemmi::CrystalSystem::Triclinic) && D > 0 && E > 0 && F > 0) { struct Flip { double kept_angle; gemmi::Mat33 basis; double d, e, f; }; const Flip flips[3] = { {uc_niggli.beta, gemmi::Mat33(-1, 0, 0, 0, 1, 0, 0, 0, -1), -D, E, -F}, @@ -444,8 +464,48 @@ std::optional SearchCharacters(const CrystalLattice &L, dou if (flip.kept_angle < 90 - NIGGLI_TYPE_BOUNDARY_DEG) continue; const auto flipped = match(L_niggli.Multiply(flip.basis), flip.d, flip.e, flip.f); - if (flipped && flipped->system != gemmi::CrystalSystem::Triclinic) - return *flipped; + if (flipped && flipped->system != gemmi::CrystalSystem::Triclinic) { + found = flipped; + break; + } + } + } + + // The characters are read on ONE reduced cell, and a lattice whose shortest vectors are degenerate + // has several: an F-centred cubic one reduces to the 60/60/60 rhombohedral cell or to a 60/60/90 + // one according to the sign of its strain, and only the first carries character 1 - the second + // reads as I-centred tetragonal or orthorhombic. Measured: of F-cubic lattices strained by 0.1 to + // 1%, 45% came back tI or oI at every strain, while the two-fold search, which does not depend on + // the basis the reduction picked, named 399 of 400 cubic F up to 0.3%. Where it finds a richer + // class than the table, take it - held to this call's own tolerances, so a search that grants no + // length equality still grants none. + // + // A tetragonal or cubic cell refused only on its length equality is, on its angles alone, the + // orthorhombic lattice of the same centring in the same basis - the class a search granting no + // length equality has to return, as the table did for the settings it could read. + if (auto lp = LePageLattice(L, std::min(angle_tolerance, LATTICE_MAX_OBLIQUITY_DEG))) { + const auto [angle, length] = ClassMetricDeparture(lp->conventional.GetUnitCell(), lp->system); + if (length > dist_tolerance && (lp->system == gemmi::CrystalSystem::Tetragonal + || lp->system == gemmi::CrystalSystem::Cubic)) + lp->system = gemmi::CrystalSystem::Orthorhombic; + if (angle <= angle_tolerance + && ClassMetricDeparture(lp->conventional.GetUnitCell(), lp->system).second <= dist_tolerance + && LatticePointGroupOrder(lp->system) > (found ? LatticePointGroupOrder(found->system) : 0) + && (!only_class || (lp->system == only_class->first && lp->centering == only_class->second))) { + int number = 0; + for (const auto &c: make_classes(D, E, F)) + if (c.system == lp->system && c.centering == lp->centering) { + number = c.number; + break; + } + return LatticeSearchResult{ + .niggli_class = number, + .primitive_reduced = lp->primitive_reduced, + .conventional = lp->conventional, + .system = lp->system, + .centering = lp->centering, + .reindex = lp->reindex, + }; } } @@ -478,6 +538,18 @@ std::optional LatticeSearchForClass(const CrystalLattice &L } +int LatticePointGroupOrder(gemmi::CrystalSystem system) { + switch (system) { + case gemmi::CrystalSystem::Monoclinic: return 2; + case gemmi::CrystalSystem::Orthorhombic: return 4; + case gemmi::CrystalSystem::Trigonal: return 6; + case gemmi::CrystalSystem::Tetragonal: return 8; + case gemmi::CrystalSystem::Hexagonal: return 12; + case gemmi::CrystalSystem::Cubic: return 24; + default: return 1; // Triclinic + } +} + bool ClassImposesLengthEquality(gemmi::CrystalSystem system) { switch (system) { case gemmi::CrystalSystem::Tetragonal: @@ -543,22 +615,7 @@ double AngleEqualityDeparture_deg(const LatticeSearchResult &sr, const CrystalLa } double MetricDeviation(const LatticeSearchResult &sr) { - const UnitCell uc = sr.conventional.GetUnitCell(); - double angle = 0.0, length = 0.0; - auto a_dev = [&](double x, double ideal) { angle = std::max(angle, std::fabs(x - ideal)); }; - auto l_dev = [&](double x, double y) { length = std::max(length, std::fabs(x - y) / (0.5 * (x + y))); }; - switch (sr.system) { - case gemmi::CrystalSystem::Monoclinic: a_dev(uc.alpha, 90); a_dev(uc.gamma, 90); break; - case gemmi::CrystalSystem::Orthorhombic: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 90); break; - case gemmi::CrystalSystem::Tetragonal: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 90); - l_dev(uc.a, uc.b); break; - case gemmi::CrystalSystem::Trigonal: - case gemmi::CrystalSystem::Hexagonal: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 120); - l_dev(uc.a, uc.b); break; - case gemmi::CrystalSystem::Cubic: a_dev(uc.alpha, 90); a_dev(uc.beta, 90); a_dev(uc.gamma, 90); - l_dev(uc.a, uc.b); l_dev(uc.b, uc.c); break; - default: break; // triclinic: nothing is asserted, nothing is violated - } + const auto [angle, length] = ClassMetricDeparture(sr.conventional.GetUnitCell(), sr.system); return std::max(angle / LATTICE_SEARCH_ANGLE_TOLERANCE_DEG, length / LATTICE_SEARCH_DIST_TOLERANCE); } diff --git a/image_analysis/lattice_search/LatticeSearch.h b/image_analysis/lattice_search/LatticeSearch.h index 7d9ef851a..5cf6588b6 100644 --- a/image_analysis/lattice_search/LatticeSearch.h +++ b/image_analysis/lattice_search/LatticeSearch.h @@ -48,6 +48,10 @@ std::optional LatticeSearchForClass(const CrystalLattice& L double dist_tolerance = LATTICE_SEARCH_DIST_TOLERANCE, double angle_tolerance_deg = LATTICE_SEARCH_ANGLE_TOLERANCE_DEG); +// Order of the lattice point group, so "lower symmetry" is a well-defined comparison +// (gemmi's enum orders Trigonal after Tetragonal, which have 6 and 8 rotations). +int LatticePointGroupOrder(gemmi::CrystalSystem system); + // Whether the conventional metric of this class asserts an equality between cell LENGTHS (a = b, and // for cubic a = b = c). Those are the promotions the relative tolerance above admits far more widely // than the spot positions resolve; every other class is carried by its angles alone. The walk run diff --git a/image_analysis/rotation_indexer/RotationIndexer.cpp b/image_analysis/rotation_indexer/RotationIndexer.cpp index 4a67037e1..6720f93d2 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.cpp +++ b/image_analysis/rotation_indexer/RotationIndexer.cpp @@ -49,20 +49,6 @@ namespace { // default iteration limit, so it only bites where the clock was biting before. constexpr int ROT_REFINE_ITERATIONS = 50; - // Order of the lattice point group, so "lower symmetry" is a well-defined comparison - // (gemmi's enum orders Trigonal after Tetragonal, which have 6 and 8 rotations). - int LatticePointGroupOrder(gemmi::CrystalSystem s) { - switch (s) { - case gemmi::CrystalSystem::Monoclinic: return 2; - case gemmi::CrystalSystem::Orthorhombic: return 4; - case gemmi::CrystalSystem::Trigonal: return 6; - case gemmi::CrystalSystem::Tetragonal: return 8; - case gemmi::CrystalSystem::Hexagonal: return 12; - case gemmi::CrystalSystem::Cubic: return 24; - default: return 1; // Triclinic - } - } - // 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. diff --git a/tests/LatticeSearchTest.cpp b/tests/LatticeSearchTest.cpp index 69e5a4ac6..72bafb701 100644 --- a/tests/LatticeSearchTest.cpp +++ b/tests/LatticeSearchTest.cpp @@ -76,6 +76,31 @@ TEST_CASE("LatticeSearch - cubic F") { CHECK(uc.gamma == Catch::Approx(90.0)); } +// An F-cubic lattice has several equally short primitive bases (60/60/60, and the 60/60/90 form that is +// also the reduced cell of an I-tetragonal lattice with c = a*sqrt(2)). Strained by a few parts per +// thousand, the reduction picks one by the sign of the strain, and the character table read on the +// 60/60/90 one named it I-centred tetragonal. +TEST_CASE("LatticeSearch - a strained cubic F lattice is cubic F in each of its primitive bases") { + const Coord a(40.0, 0.0, 0.0), b(0.05, 40.08, 0.0), c(0.0, -0.04, 39.94); + const CrystalLattice bases[] = { + CrystalLattice((b + c) * 0.5, (a + c) * 0.5, (a + b) * 0.5), + CrystalLattice((a + b) * 0.5, (a - b) * 0.5, (a + c) * 0.5), + CrystalLattice((a + b) * 0.5, (b - a) * 0.5, (a + c) * 0.5), + }; + for (const auto &L : bases) { + const auto res = LatticeSearch(L); + CHECK(res.system == gemmi::CrystalSystem::Cubic); + CHECK(res.centering == 'F'); + const auto uc = res.conventional.GetUnitCell(); + CHECK(uc.a == Catch::Approx(40.0).epsilon(0.01)); + CHECK(uc.b == Catch::Approx(40.0).epsilon(0.01)); + CHECK(uc.c == Catch::Approx(40.0).epsilon(0.01)); + CHECK(uc.alpha == Catch::Approx(90.0).margin(0.5)); + CHECK(uc.beta == Catch::Approx(90.0).margin(0.5)); + CHECK(uc.gamma == Catch::Approx(90.0).margin(0.5)); + } +} + TEST_CASE("LatticeSearch - cubic P") { // Simple cubic P, a=30 const double a = 30.0;