diff --git a/image_analysis/lattice_search/LatticeSearch.cpp b/image_analysis/lattice_search/LatticeSearch.cpp index 1f3ae4f7..fbeb22a6 100644 --- a/image_analysis/lattice_search/LatticeSearch.cpp +++ b/image_analysis/lattice_search/LatticeSearch.cpp @@ -6,6 +6,7 @@ #include #include +#include struct NiggliClass { int number; @@ -39,7 +40,9 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance double E = g_vec.eta / 2; double F = g_vec.zeta / 2; - std::vector niggli_classes = { + // D, E, F are parameters so the table can also be built for the type-flipped setting below. + auto make_classes = [&](double D, double E, double F) { + return std::vector{ { 1, 1, true, true, A / 2, A / 2, A / 2, false, false, @@ -301,48 +304,69 @@ LatticeSearchResult LatticeSearch(const CrystalLattice &L, double dist_tolerance {1, 0, 0, 0, 1, 0, 0, 0, 1}, gemmi::CrystalSystem::Triclinic, 'P' } + }; }; - const auto uc_reduced = L_niggli.GetUnitCell(); + auto match = [&](const CrystalLattice &latt, double D, double E, double F) + -> std::optional { + const auto uc_reduced = latt.GetUnitCell(); + for (const auto &c: make_classes(D, E, F)) { + if (c.type == 1 && uc_reduced.beta >= 90 - angle_tolerance ) + continue; - for (const auto &c: niggli_classes) { - if (c.type == 1 && uc_reduced.beta >= 90 - angle_tolerance ) - continue; + bool ok = true; - bool ok = true; + if (c.cond_AB && fabs((uc_reduced.a - uc_reduced.b) / (0.5 * (uc_reduced.a + uc_reduced.b))) > dist_tolerance) + ok = false; + if (c.cond_BC && fabs((uc_reduced.b - uc_reduced.c) / (0.5 * (uc_reduced.b + uc_reduced.c))) > dist_tolerance) + ok = false; - if (c.cond_AB && fabs((uc_reduced.a - uc_reduced.b) / (0.5 * (uc_reduced.a + uc_reduced.b))) > dist_tolerance) - ok = false; - if (c.cond_BC && fabs((uc_reduced.b - uc_reduced.c) / (0.5 * (uc_reduced.b + uc_reduced.c))) > dist_tolerance) - ok = false; + double expected_alpha = acos(c.cond_D / sqrt(B*C)) * 180 / PI; + double expected_beta = acos(c.cond_E / sqrt(A*C)) * 180 / PI; + double expected_gamma = acos(c.cond_F / sqrt(A*B)) * 180 / PI; - double expected_alpha = acos(c.cond_D / sqrt(B*C)) * 180 / PI; - double expected_beta = acos(c.cond_E / sqrt(A*C)) * 180 / PI; - double expected_gamma = acos(c.cond_F / sqrt(A*B)) * 180 / PI; + if (fabs(expected_alpha - uc_reduced.alpha) > angle_tolerance) + ok = false; + if (fabs(expected_beta - uc_reduced.beta) > angle_tolerance) + ok = false; + if (fabs(expected_gamma - uc_reduced.gamma) > angle_tolerance) + ok = false; - if (fabs(expected_alpha - uc_reduced.alpha) > angle_tolerance) - ok = false; - if (fabs(expected_beta - uc_reduced.beta) > angle_tolerance) - ok = false; - if (fabs(expected_gamma - uc_reduced.gamma) > angle_tolerance) - ok = false; + double tmp1 = 2.0 * fabs(D + E + F); + double tmp2 = A + B; - double tmp1 = 2.0 * fabs(D + E + F); - double tmp2 = A + B; + if (c.cond_DEF && fabs((tmp1 - tmp2) / (0.5 * (tmp1 + tmp2))) > dist_tolerance) + ok = false; - if (c.cond_DEF && fabs((tmp1 - tmp2) / (0.5 * (tmp1 + tmp2))) > dist_tolerance) - ok = false; - - if (ok) { - return LatticeSearchResult{ - .niggli_class = c.number, - .primitive_reduced = L_niggli, - .conventional = L_niggli.Multiply(c.reindex), - .system = c.system, - .centering = c.centering, - .reindex = c.reindex, - }; + if (ok) { + return LatticeSearchResult{ + .niggli_class = c.number, + .primitive_reduced = latt, + .conventional = latt.Multiply(c.reindex), + .system = c.system, + .centering = c.centering, + .reindex = c.reindex, + }; + } } + return std::nullopt; + }; + + if (auto found = match(L_niggli, D, E, F)) + return *found; + + // A reduced cell whose beta is 90 to within the angle tolerance 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 the + // last digits of whatever refinement produced it. The type-1 characters are skipped for such a cell + // (just above) and the type-2 ones are stated for the obtuse setting, so an acute cell can match none + // of them and comes back triclinic. Present it in the obtuse setting and try once more - negating a + // and c keeps the lattice and beta and turns alpha and gamma into their supplements. Measured on a + // C-centred monoclinic crystal whose reduced beta sits 0.07 deg from 90: its centring was read or + // missed according to the sign of that 0.07 deg, and with it the space group of the whole run. + if (L_niggli.GetUnitCell().beta >= 90 - angle_tolerance && D > 0 && E > 0 && F > 0) { + const gemmi::Mat33 obtuse(-1, 0, 0, 0, 1, 0, 0, 0, -1); + if (auto found = match(L_niggli.Multiply(obtuse), -D, E, -F)) + return *found; } return LatticeSearchResult{