mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 67 Wiki Activity

Lattice search: a reduced cell sitting on the Niggli type boundary keeps its centring

Browse Source 
The 44 lattice characters come in two types - all-acute reduced cells and all-obtuse ones - and
the search skips the type-1 characters whenever the reduced beta is within the angle tolerance of
90 degrees, because at that point the two types are no longer distinguishable and only the type-2
statement of a character is safe to test. But a cell whose reduced beta really is 90 within
tolerance can itself reduce EITHER way, and an all-acute one then matches no character at all: the
type-1 ones were skipped and the type-2 ones are written for the obtuse setting. It falls through
to triclinic, and a centred lattice loses its centring.

Measured on a C-centred monoclinic crystal whose reduced beta sits 0.07 degrees from 90. Its
free-refined cell came back acute or obtuse depending on the last digits of the refinement - a
sub-pixel change in the beam centre was enough - and with it the lattice was read as C-centred
monoclinic or as triclinic, and the run merged in C2 or in P1 with twice the asymmetric unit. In
the obtuse setting the character matches with residuals of 0.01 to 0.14 degrees against a 3 degree
tolerance, at every geometry tried; there was never any doubt about the lattice, only about which
side of the boundary the reduction landed on.

So when nothing matched and the cell is acute with beta at the boundary, present it in the obtuse
setting and match once more. Negating a and c leaves the lattice and beta alone and turns alpha
and gamma into their supplements. The second attempt runs only where the first found nothing, so
no lattice that is classified today can be re-classified by this.
This commit is contained in:
leonarski_f
2026-08-12 17:02:36 +02:00
parent 8874a788e6
commit 996d8e7d41
1 changed files with 57 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/lattice_search/LatticeSearch.cpp
+57 -33
View File
@@ -6,6 +6,7 @@
#include <gemmi/cellred.hpp>
#include <cmath>
#include <optional>
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<NiggliClass> 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<NiggliClass>{
{
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<LatticeSearchResult> {
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{