CrystalLattice: Handle Monoclinic cell to enforce obtuse beta

This commit is contained in:
2026-03-03 14:36:08 +01:00
parent f3e0a15d26
commit 2e08bff4df
3 changed files with 6 additions and 34 deletions

View File

@@ -139,35 +139,6 @@ void CrystalLattice::ReorderABEqual() {
}
void CrystalLattice::ReorderMonoclinic() {
float alpha = angle_deg(vec[1], vec[2]);
float beta = angle_deg(vec[0], vec[2]);
float gamma = angle_deg(vec[0], vec[1]);
float da = std::abs(alpha - 90.0f);
float db = std::abs(beta - 90.0f);
float dg = std::abs(gamma - 90.0f);
Coord a = vec[0];
Coord b = vec[1];
Coord c = vec[2];
if (da > db && da > dg) {
// alpha most different -> [b, c, a]
vec[0] = b;
vec[1] = a;
vec[2] = c;
} else if (dg > db && dg > da) {
// gamma most different -> [c, a, b]
vec[0] = a;
vec[1] = c;
vec[2] = b;
} // else beta most different -> keep [a, b, c]
if (vec[0].Length() > vec[2].Length())
std::swap(vec[0], vec[2]);
FixHandeness();
// Enforce obtuse beta (>= 90°). Beta is the angle between a and c.
// Flip signs of a and b simultaneously to keep handedness and lengths unchanged,
// which maps beta -> 180° - beta.

View File

@@ -82,6 +82,9 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetry(Data
.crystal_system = sym_result.system
};
outcome.lattice_candidate = sym_result.conventional;
if (sym_result.system == gemmi::CrystalSystem::Monoclinic)
outcome.lattice_candidate->ReorderMonoclinic();
}
return outcome;

View File

@@ -67,10 +67,8 @@ TEST_CASE("CrystalLattice_Sort") {
TEST_CASE("CrystalLattice_ReorderMonoclinic") {
std::vector<CrystalLattice> latt = {
{85,70,60, 85, 90, 90},
{60,70,85, 90, 90, 85},
{60,85,70, 90, 85, 90},
{70,60,85, 90, 90, 85}
{60, 85, 70, 90, 70, 90},
{60, 85, 70, 90, 110, 90},
};
for (const auto &l_in :latt) {
CrystalLattice l = l_in;
@@ -79,7 +77,7 @@ TEST_CASE("CrystalLattice_ReorderMonoclinic") {
CHECK(l.Vec0().Length() == Catch::Approx(60));
CHECK(l.Vec1().Length() == Catch::Approx(85));
CHECK(l.Vec2().Length() == Catch::Approx(70));
CHECK(l.GetUnitCell().beta == Catch::Approx(95.0));
CHECK(l.GetUnitCell().beta == Catch::Approx(110.0));
}
}