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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EFEJG6WBQv8th4UJFNe53N
(cherry picked from commit 6ca00e927d664e870515c04164defa81d8a18725)
This commit is contained in:
2026-08-31 07:17:00 +02:00
parent 527a5187f4
commit e0ff51e2a0
3 changed files with 32 additions and 2 deletions
+22
View File
@@ -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));
}