From e0ff51e2a0b58d64432b6055e72f587d1c7bb034 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 30 Aug 2026 17:19:11 +0200 Subject: [PATCH] lattice: transpose the change of basis to a primitive cell, as everywhere else gemmi states centred_to_primitive as an operator on COORDINATES and CrystalLattice::Multiply combines BASIS VECTORS, so the matrix has to be transposed on the way in - as it already is at both of the other places a gemmi Op::Rot reaches Multiply, one of them in this same file. A, B, C, I and F are symmetric matrices, so for them the transpose is a no-op and the omission never showed. R and H are not. Measured on an R-centred hexagonal lattice, ToPrimitive('R') returned 59.5 81.7 43.3 / 145.6 124.5 46.7 where the rhombohedral primitive cell is 49.3 49.3 49.3 / 60.9 60.9 60.9. What hid it is that a determinant is unchanged by transposition, so the VOLUME came out right - and most callers only take the volume. It is not only cosmetic: the result feeds the re-seating path that puts a lattice into a space group the user fixed by hand, so an R-centred lattice was handed the classifier a "primitive" cell that is not that lattice - broken for exactly the centring whose setting most needs re-seating. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N (cherry picked from commit 6ca00e927d664e870515c04164defa81d8a18725) --- common/CrystalLattice.cpp | 11 +++++++++-- docs/CHANGELOG.md | 1 + tests/LatticeSearchTest.cpp | 22 ++++++++++++++++++++++ 3 files changed, 32 insertions(+), 2 deletions(-) diff --git a/common/CrystalLattice.cpp b/common/CrystalLattice.cpp index 7a79630ae..b5c225cde 100644 --- a/common/CrystalLattice.cpp +++ b/common/CrystalLattice.cpp @@ -200,14 +200,21 @@ CrystalLattice CrystalLattice::FromPrimitive(char centering) const { if (centering == 'P') return *this; - return Multiply(gemmi::rot_as_mat33(gemmi::centred_to_primitive(centering)).inverse()); + // Transposed for the same reason as ToPrimitive below. + return Multiply(gemmi::rot_as_mat33(gemmi::centred_to_primitive(centering)).transpose().inverse()); } CrystalLattice CrystalLattice::ToPrimitive(char centering) const { if (centering == 'P') return *this; - return Multiply(gemmi::rot_as_mat33(gemmi::centred_to_primitive(centering))); + // gemmi states the change of basis as an operator on COORDINATES, and Multiply combines BASIS + // VECTORS, so the matrix has to be transposed - as it already is everywhere else a gemmi Op::Rot + // reaches Multiply. A, B, C, I and F are symmetric, so the transpose is a no-op for them and the + // omission never showed; R and H are not, and without it an R-centred lattice was handed back a + // "primitive" cell that is not that lattice. Its VOLUME was right either way, which is what hid + // this: a determinant does not change under transposition, and most callers only take the volume. + return Multiply(gemmi::rot_as_mat33(gemmi::centred_to_primitive(centering)).transpose()); } void CrystalLattice::Regularize(const gemmi::CrystalSystem &input) { diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 557ae0886..51799229f 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -9,6 +9,7 @@ * `rugnux` writes the unmerged MTZ `_unmerged.mtz` by default; `--no-export-unmerged` skips it. * Every rotation run that determines its own space group also writes `_P1.mtz`, the same observations merged in P1, so a wrong space group can be re-merged, re-solved or re-refined without processing the images again; `--no-p1-crosscheck` declines it, and a run given `-S` writes nothing because its space group's centring absences were never integrated. * The space-group search names the setting the data support, so a screw or a 2-fold on the a or c axis is reported as such instead of costing the crystal its space group. +* The primitive cell of an R- or H-centred lattice is that lattice: the change of basis was applied without the transposition the convention needs, which left the volume right and the cell wrong, and reached the path that re-seats a lattice into a space group fixed by hand. * A centred monoclinic lattice that reduces to the mI form keeps its centring: the character naming that reduced form was missing from the table, so such a cell was reported as triclinic. * A lattice keeps the symmetry it has whatever orientation it was refined in: the cell reduction now judges a structurally-zero scalar product against the size of the cell rather than against a fixed tolerance eight decades smaller, so a body-centred tetragonal lattice is no longer read as C-centred monoclinic on most rotations of itself. * A run whose cell metric admits more rotational symmetry than the space group it adopted says so, naming both, so a symmetry the intensities were too weak to confirm is visible rather than silent. diff --git a/tests/LatticeSearchTest.cpp b/tests/LatticeSearchTest.cpp index 216e8201e..bd2fcecfc 100644 --- a/tests/LatticeSearchTest.cpp +++ b/tests/LatticeSearchTest.cpp @@ -574,3 +574,25 @@ TEST_CASE("LatticeSearch - a centred monoclinic lattice that reduces to the mI f CHECK(std::fabs(res.conventional.CalcVolume()) == Catch::Approx(2 * std::fabs(L.CalcVolume())).epsilon(1e-4)); } + +// The change of basis to a primitive cell is stated by gemmi as an operator on COORDINATES, while +// CrystalLattice::Multiply combines BASIS VECTORS, so it has to be transposed. A, B, C, I and F are +// symmetric and never showed the omission; R and H are not. An R-centred lattice is the case that +// matters, because it is the centring whose setting most often has to be re-seated. +TEST_CASE("CrystalLattice::ToPrimitive gives an R-centred lattice its rhombohedral primitive cell") { + const double a = 50.0, c = 120.0; + const CrystalLattice hex(a, a, c, 90, 90, 120); + const auto prim = hex.ToPrimitive('R').GetUnitCell(); + // A rhombohedral primitive cell: three equal edges, three equal angles, a third of the volume. + CHECK(prim.a == Catch::Approx(prim.b).epsilon(1e-5)); + CHECK(prim.b == Catch::Approx(prim.c).epsilon(1e-5)); + CHECK(prim.alpha == Catch::Approx(prim.beta).epsilon(1e-5)); + CHECK(prim.beta == Catch::Approx(prim.gamma).epsilon(1e-5)); + CHECK(std::fabs(hex.ToPrimitive('R').CalcVolume()) + == Catch::Approx(std::fabs(hex.CalcVolume()) / 3.0).epsilon(1e-4)); + // ...and it goes back to the hexagonal cell it came from. + const auto back = hex.ToPrimitive('R').FromPrimitive('R').GetUnitCell(); + CHECK(back.a == Catch::Approx(a).epsilon(1e-4)); + CHECK(back.c == Catch::Approx(c).epsilon(1e-4)); + CHECK(back.gamma == Catch::Approx(120.0).epsilon(1e-4)); +}