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;