leonarski_fandClaude Opus 5
b877d61a42
The departure from a = b that decides whether the metric-symmetry arms run at all was read by pushing the free (triclinic) refinement through the search result's `reindex`. That matrix is stale for a whole population: where the walk matches a metrically hexagonal lattice on its C-centred orthorhombic character, RotationIndexer re-expresses `conventional` in hexagonal axes but leaves `reindex` describing the setting it replaced. The C-centred setting has b = a*sqrt(3) BY CONSTRUCTION, so a against b there is 2(sqrt(3)-1)/(sqrt(3)+1) = 53.6 % apart on every hexagonal lattice there is, whatever the crystal. So the gate fired on the whole hexagonal and trigonal population - eighteen sets in one corpus run, all reporting 53.1-53.7 % - and the arms then decided two classes that fit equally well on differences in the third or fourth significant figure. Three of those coin flips landed on the ortho-hexagonal supercell and took a P 61 2 2, a P 65 2 2 and a P 31 with them. The departure is now computed by LengthEqualityDeparture, which derives the primitive-to-conventional map from the pair of cells the search result carries instead of trusting `reindex`, and reads a against b in that basis. The two cells describe one lattice, so the map is the integer relabelling of its basis vectors. The measurement no longer depends on which character the walk happened to match. Re-measured, the hexagonal and trigonal sets report 0.00-0.41 %, the band the tetragonal and cubic proteins were already in, and the gate does not fire on any of them; the genuine pseudo-symmetric small-molecule cases stand alone at 0.98, 1.85 and 2.61 %. The three lost groups come back, and the sets that merely tolerated the noise are bit-identical. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
711 lines
30 KiB
C++
711 lines
30 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#include <catch2/catch_all.hpp>
|
|
#include "../common/CrystalLattice.h"
|
|
#include "../common/Coord.h"
|
|
#include "../common/UnitCell.h"
|
|
#include "../image_analysis/lattice_search/LatticeSearch.h"
|
|
#include "gemmi/symmetry.hpp"
|
|
#include <cmath>
|
|
|
|
// Helper: check near-equality of unit cell parameters
|
|
static void check_uc(const UnitCell& uc, double a, double b, double c,
|
|
double alpha, double beta, double gamma,
|
|
double eps_len = 1e-6, double eps_ang = 1e-4) {
|
|
CHECK(uc.a == Catch::Approx(a).margin(eps_len));
|
|
CHECK(uc.b == Catch::Approx(b).margin(eps_len));
|
|
CHECK(uc.c == Catch::Approx(c).margin(eps_len));
|
|
CHECK(uc.alpha == Catch::Approx(alpha).margin(eps_ang));
|
|
CHECK(uc.beta == Catch::Approx(beta ).margin(eps_ang));
|
|
CHECK(uc.gamma == Catch::Approx(gamma).margin(eps_ang));
|
|
}
|
|
|
|
|
|
TEST_CASE("LatticeSearch - cubic I") {
|
|
// Build a body-centered cubic cell with a=40:
|
|
// primitive basis vectors (conventional I cubic primitive):
|
|
// p1 = (0, a/2, a/2), p2 = (a/2, 0, a/2), p3 = (a/2, a/2, 0)
|
|
const double a = 40.0;
|
|
CrystalLattice L(
|
|
Coord(a, 0, 0),
|
|
Coord(0, a, 0),
|
|
Coord(0, 0, a)
|
|
);
|
|
L = L.ToPrimitive('I');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Cubic);
|
|
CHECK(res.centering == 'I');
|
|
|
|
// Conventional cubic I should have equal edges and 90° angles
|
|
auto uc = res.conventional.GetUnitCell();
|
|
CHECK(uc.a == Catch::Approx( a )); // In this construction, conventional a matches given a
|
|
CHECK(uc.b == Catch::Approx( a ));
|
|
CHECK(uc.c == Catch::Approx( a ));
|
|
CHECK(uc.alpha == Catch::Approx(90.0));
|
|
CHECK(uc.beta == Catch::Approx(90.0));
|
|
CHECK(uc.gamma == Catch::Approx(90.0));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - cubic F") {
|
|
// Build a body-centered cubic cell with a=40:
|
|
// primitive basis vectors (conventional I cubic primitive):
|
|
// p1 = (0, a/2, a/2), p2 = (a/2, 0, a/2), p3 = (a/2, a/2, 0)
|
|
const double a = 40.0;
|
|
CrystalLattice L(
|
|
Coord(a, 0, 0),
|
|
Coord(0, a, 0),
|
|
Coord(0, 0, a)
|
|
);
|
|
L = L.ToPrimitive('F');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Cubic);
|
|
CHECK(res.centering == 'F');
|
|
|
|
// Conventional cubic I should have equal edges and 90° angles
|
|
auto uc = res.conventional.GetUnitCell();
|
|
CHECK(uc.a == Catch::Approx( a )); // In this construction, conventional a matches given a
|
|
CHECK(uc.b == Catch::Approx( a ));
|
|
CHECK(uc.c == Catch::Approx( a ));
|
|
CHECK(uc.alpha == Catch::Approx(90.0));
|
|
CHECK(uc.beta == Catch::Approx(90.0));
|
|
CHECK(uc.gamma == Catch::Approx(90.0));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - cubic P") {
|
|
// Simple cubic P, a=30
|
|
const double a = 30.0;
|
|
CrystalLattice L(
|
|
Coord(a,0,0),
|
|
Coord(0,a,0),
|
|
Coord(0,0,a)
|
|
);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Cubic);
|
|
CHECK(res.centering == 'P');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, a, a, a, 90.0, 90.0, 90.0, 1e-6, 1e-4);
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - tetragonal I") {
|
|
// Build a body-centered cubic cell with a=40:
|
|
// primitive basis vectors (conventional I cubic primitive):
|
|
// p1 = (0, a/2, a/2), p2 = (a/2, 0, a/2), p3 = (a/2, a/2, 0)
|
|
const double a = 40.0;
|
|
const double b = 34.0;
|
|
CrystalLattice L(
|
|
Coord(a, 0, 0),
|
|
Coord(0, a, 0),
|
|
Coord(0, 0, b)
|
|
);
|
|
L = L.ToPrimitive('I');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Tetragonal);
|
|
CHECK(res.centering == 'I');
|
|
|
|
// Conventional cubic I should have equal edges and 90° angles
|
|
auto uc = res.conventional.GetUnitCell();
|
|
CHECK(uc.a == Catch::Approx( a )); // In this construction, conventional a matches given a
|
|
CHECK(uc.b == Catch::Approx( a ));
|
|
CHECK(uc.c == Catch::Approx( b ));
|
|
CHECK(uc.alpha == Catch::Approx(90.0));
|
|
CHECK(uc.beta == Catch::Approx(90.0));
|
|
CHECK(uc.gamma == Catch::Approx(90.0));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - tetragonal I - v2") {
|
|
// Build a body-centered cubic cell with a=40:
|
|
// primitive basis vectors (conventional I cubic primitive):
|
|
// p1 = (0, a/2, a/2), p2 = (a/2, 0, a/2), p3 = (a/2, a/2, 0)
|
|
const double a = 40.0;
|
|
const double b = 54.0;
|
|
CrystalLattice L(
|
|
Coord(a, 0, 0),
|
|
Coord(0, a, 0),
|
|
Coord(0, 0, b)
|
|
);
|
|
L = L.ToPrimitive('I');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Tetragonal);
|
|
CHECK(res.centering == 'I');
|
|
|
|
// Conventional cubic I should have equal edges and 90° angles
|
|
auto uc = res.conventional.GetUnitCell();
|
|
CHECK(uc.a == Catch::Approx( a )); // In this construction, conventional a matches given a
|
|
CHECK(uc.b == Catch::Approx( a ));
|
|
CHECK(uc.c == Catch::Approx( b ));
|
|
CHECK(uc.alpha == Catch::Approx(90.0));
|
|
CHECK(uc.beta == Catch::Approx(90.0));
|
|
CHECK(uc.gamma == Catch::Approx(90.0));
|
|
}
|
|
|
|
// Tetragonal P: a=b!=c, all angles 90, P-centering
|
|
TEST_CASE("LatticeSearch - tetragonal P") {
|
|
const double a = 37.0, c = 59.0;
|
|
CrystalLattice L(
|
|
Coord(a,0,0),
|
|
Coord(0,a,0),
|
|
Coord(0,0,c)
|
|
);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Tetragonal);
|
|
CHECK(res.centering == 'P');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, a, a, c, 90.0, 90.0, 90.0, 1e-2, 1e-2);
|
|
}
|
|
|
|
// Orthorhombic F: all angles 90, unequal edges, F-centering
|
|
TEST_CASE("LatticeSearch - orthorhombic F") {
|
|
const double a = 35.0, b = 41.0, c = 57.0;
|
|
CrystalLattice conv(a,b,c, 90.0,90.0,90.0);
|
|
|
|
CrystalLattice L = conv.ToPrimitive('F');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'F');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-1, 1e-2);
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - orthorhombic F - permutation 1") {
|
|
const double a = 41.0, b = 57.0, c = 35.0;
|
|
CrystalLattice conv(a,b,c, 90.0,90.0,90.0);
|
|
|
|
CrystalLattice L = conv.ToPrimitive('F');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'F');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, c, a, b, 90.0, 90.0, 90.0, 1e-1, 1e-2);
|
|
}
|
|
|
|
// Orthorhombic C: all angles 90, unequal edges, C-centering
|
|
TEST_CASE("LatticeSearch - orthorhombic C") {
|
|
const double a = 35.0, b = 41.0, c = 57.0;
|
|
CrystalLattice conv(a,b,c, 90.0,90.0,90.0);
|
|
CrystalLattice L = conv.ToPrimitive('C');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'C');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-1, 1e-2);
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - orthorhombic I") {
|
|
const double a = 35.0, b = 41.0, c = 57.0;
|
|
CrystalLattice conv(a,b,c, 90.0,90.0,90.0);
|
|
CrystalLattice L = conv.ToPrimitive('I');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'I');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-2, 1e-2);
|
|
}
|
|
|
|
|
|
TEST_CASE("LatticeSearch - orthorhombic I - permutation1") {
|
|
const double a = 57.0, b = 41.0, c = 35.0;
|
|
CrystalLattice conv(a,b,c, 90.0,90.0,90.0);
|
|
CrystalLattice L = conv.ToPrimitive('I');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'I');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, c, b, a, 90.0, 90.0, 90.0, 1e-2, 1e-2);
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - orthorhombic I - permutation2") {
|
|
const double a = 41.0, b = 57.0, c = 35.0;
|
|
CrystalLattice conv(a,b,c, 90.0,90.0,90.0);
|
|
CrystalLattice L = conv.ToPrimitive('I');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'I');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, c, a, b, 90.0, 90.0, 90.0, 1e-2, 1e-2);
|
|
}
|
|
|
|
// A character states its scalar products as fractions of A, B and C, and the three C-centred
|
|
// monoclinic ones (28, 29, 30) state one of them as 2*D or 2*E - twice a cosine. The cosine that
|
|
// implies leaves [-1,1] as soon as the cell's own angle is far enough from 90, and the character is
|
|
// then geometrically impossible for that metric. This cell is triclinic; character 28 asks it for a
|
|
// gamma whose cosine is 1.127.
|
|
TEST_CASE("LatticeSearch - an impossible character is not a match") {
|
|
CrystalLattice L(30.0, 35.0, 40.0, 65.0, 70.0, 70.0);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Triclinic);
|
|
}
|
|
|
|
// An exact I-centred orthorhombic lattice whose reduced cell comes out all-acute with gamma at 90 -
|
|
// ON the boundary between the two Niggli types, where the reduction may present either. Character 42
|
|
// is stated for the obtuse setting, and only the flip that keeps gamma reaches it. Both defects have
|
|
// to be gone: without the impossible-character fix this metric matches character 28 and never gets
|
|
// as far as the retry, and without the gamma flip the retry does not have the setting it needs.
|
|
TEST_CASE("LatticeSearch - orthorhombic I on the type boundary in gamma") {
|
|
const double a = 45.0, b = 50.0, c = 80.0;
|
|
CrystalLattice conv(a, b, c, 90.0, 90.0, 90.0);
|
|
CrystalLattice L = conv.ToPrimitive('I');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'I');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-2, 1e-2);
|
|
}
|
|
|
|
// Orthorhombic P: all angles 90, unequal edges, P-centering
|
|
TEST_CASE("LatticeSearch - orthorhombic P") {
|
|
const double a = 35.0, b = 41.0, c = 57.0;
|
|
CrystalLattice L(a,b,c, 90.0,90.0,90.0);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(res.centering == 'P');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-6, 1e-4);
|
|
}
|
|
|
|
// Hexagonal P: a=b!=c, alpha=beta=90, gamma=120, P-centering
|
|
TEST_CASE("LatticeSearch - hexagonal P") {
|
|
const double a = 30.0, c = 48.0;
|
|
CrystalLattice L(
|
|
Coord(a, 0, 0),
|
|
Coord(-a/2, a*std::sqrt(3)/2, 0),
|
|
Coord(0, 0, c)
|
|
);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Hexagonal);
|
|
CHECK(res.centering == 'P');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
CHECK(uc.a == Catch::Approx(a).margin(1e-2));
|
|
CHECK(uc.b == Catch::Approx(a).margin(1e-2));
|
|
CHECK(uc.c == Catch::Approx(c).margin(1e-2));
|
|
CHECK(uc.alpha == Catch::Approx(90.0).margin(1e-2));
|
|
CHECK(uc.beta == Catch::Approx(90.0).margin(1e-2));
|
|
CHECK(uc.gamma == Catch::Approx(120.0).margin(1e-2));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - monoclinic C (unique b)") {
|
|
const double a = 50.0, b = 60.0, c = 70.0;
|
|
const double alpha = 90.0, beta = 96.0, gamma = 90.0;
|
|
CrystalLattice conv(a,b,c, alpha,beta,gamma);
|
|
auto L = conv.ToPrimitive('C');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Monoclinic);
|
|
CHECK(res.centering == 'C');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
// Check right angles at alpha,gamma and non-90 beta; lengths comparable
|
|
CHECK(std::fabs(uc.alpha - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.gamma - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.beta - beta) < 1e-2);
|
|
// Lengths should match within small tolerance
|
|
CHECK(uc.a == Catch::Approx(a).margin(1e-2));
|
|
CHECK(uc.b == Catch::Approx(b).margin(1e-2));
|
|
CHECK(uc.c == Catch::Approx(c).margin(1e-2));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - monoclinic C (unique b) - v2") {
|
|
const double a = 71.0, b = 35.0, c = 90.0;
|
|
const double alpha = 90.0, beta = 96.0, gamma = 90.0;
|
|
CrystalLattice conv(a,b,c, alpha,beta,gamma);
|
|
auto L = conv.ToPrimitive('C');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Monoclinic);
|
|
CHECK(res.centering == 'C');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
// Check right angles at alpha,gamma and non-90 beta; lengths comparable
|
|
CHECK(std::fabs(uc.alpha - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.gamma - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.beta - beta) < 1e-2);
|
|
// Lengths should match within small tolerance
|
|
CHECK(uc.a == Catch::Approx(a).margin(1e-2));
|
|
CHECK(uc.b == Catch::Approx(b).margin(1e-2));
|
|
CHECK(uc.c == Catch::Approx(c).margin(1e-2));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - monoclinic C (unique a)") {
|
|
const double a = 60.0, b = 50.0, c = 70.0;
|
|
const double alpha = 96.0, beta = 90.0, gamma = 90.0;
|
|
CrystalLattice conv(a,b,c, alpha,beta,gamma);
|
|
auto L = conv.ToPrimitive('C');
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Monoclinic);
|
|
CHECK(res.centering == 'C');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
// Check right angles at alpha,gamma and non-90 beta; lengths comparable
|
|
CHECK(std::fabs(uc.alpha - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.gamma - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.beta - alpha) < 1e-2);
|
|
// Lengths should match within small tolerance
|
|
CHECK(uc.a == Catch::Approx(b).margin(1e-2));
|
|
CHECK(uc.b == Catch::Approx(a).margin(1e-2));
|
|
CHECK(uc.c == Catch::Approx(c).margin(1e-2));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - monoclinic P (unique b)") {
|
|
const double a = 50.0, b = 60.0, c = 70.0;
|
|
const double alpha = 90.0, beta = 96.0, gamma = 90.0;
|
|
CrystalLattice conv(a,b,c, alpha,beta,gamma);
|
|
|
|
auto res = LatticeSearch(conv, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Monoclinic);
|
|
CHECK(res.centering == 'P');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
// Check right angles at alpha,gamma and non-90 beta; lengths comparable
|
|
CHECK(std::fabs(uc.alpha - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.gamma - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.beta - beta) < 1e-2);
|
|
// Lengths should match within small tolerance
|
|
CHECK(uc.a == Catch::Approx(a).margin(1e-2));
|
|
CHECK(uc.b == Catch::Approx(b).margin(1e-2));
|
|
CHECK(uc.c == Catch::Approx(c).margin(1e-2));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - monoclinic P (unique b) - v2") {
|
|
const double a = 90.0, b = 35.0, c = 71.0;
|
|
const double alpha = 90.0, beta = 96.0, gamma = 90.0;
|
|
CrystalLattice conv(a,b,c, alpha,beta,gamma);
|
|
|
|
auto res = LatticeSearch(conv, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Monoclinic);
|
|
CHECK(res.centering == 'P');
|
|
|
|
auto uc = res.conventional.GetUnitCell();
|
|
// Check right angles at alpha,gamma and non-90 beta; lengths comparable
|
|
CHECK(std::fabs(uc.alpha - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.gamma - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.beta - beta) < 1e-2);
|
|
// Lengths should match within small tolerance
|
|
CHECK(uc.a == Catch::Approx(c).margin(1e-2));
|
|
CHECK(uc.b == Catch::Approx(b).margin(1e-2));
|
|
CHECK(uc.c == Catch::Approx(a).margin(1e-2));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - triclinic P") {
|
|
// General triclinic primitive cell
|
|
CrystalLattice L(33.1, 41.7, 52.3, 89.1, 85.0, 76.3);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
// System should be triclinic, centering P, and conventional equals some standardized primitive
|
|
CHECK(res.system == gemmi::CrystalSystem::Triclinic);
|
|
CHECK(res.centering == 'P');
|
|
|
|
// The conventional cell should be metric-equivalent to input. We verify only the system and centering here.
|
|
// Reduced primitive must be non-singular
|
|
auto uc_red = res.primitive_reduced.GetUnitCell();
|
|
CHECK(uc_red.a > 0);
|
|
CHECK(uc_red.b > 0);
|
|
CHECK(uc_red.c > 0);
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - triclinic P - v2") {
|
|
// General triclinic primitive cell
|
|
CrystalLattice L(33.1, 41.7, 52.3, 100, 92, 115);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
// System should be triclinic, centering P, and conventional equals some standardized primitive
|
|
CHECK(res.system == gemmi::CrystalSystem::Triclinic);
|
|
CHECK(res.centering == 'P');
|
|
|
|
// The conventional cell should be metric-equivalent to input. We verify only the system and centering here.
|
|
// Reduced primitive must be non-singular
|
|
auto uc_red = res.primitive_reduced.GetUnitCell();
|
|
CHECK(uc_red.a > 0);
|
|
CHECK(uc_red.b > 0);
|
|
CHECK(uc_red.c > 0);
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - trigonal R") {
|
|
const double a = 32.0;
|
|
const double alpha = 80.0;
|
|
|
|
// Build rhombohedral in rhombohedral setting (primitive axes a=b=c, alpha=beta=gamma)
|
|
CrystalLattice L(a, a, a, alpha, alpha, alpha);
|
|
|
|
auto res = LatticeSearch(L, 1e-6);
|
|
|
|
CHECK(res.system == gemmi::CrystalSystem::Trigonal);
|
|
CHECK(res.centering == 'R');
|
|
|
|
auto uc_red = res.conventional.GetUnitCell();
|
|
CHECK(uc_red.alpha == Catch::Approx(90).margin(1e-2));
|
|
CHECK(uc_red.beta == Catch::Approx(90).margin(1e-2));
|
|
CHECK(uc_red.gamma == Catch::Approx(120).margin(1e-2));
|
|
|
|
auto uc_prim = res.primitive_reduced.GetUnitCell();
|
|
CHECK(uc_prim.alpha == Catch::Approx(alpha).margin(1e-2));
|
|
CHECK(uc_prim.beta == Catch::Approx(alpha).margin(1e-2));
|
|
CHECK(uc_prim.gamma == Catch::Approx(alpha).margin(1e-2));
|
|
}
|
|
|
|
// The class-filtered walk: the same table, restricted to one Bravais class. A tetragonal-P lattice is
|
|
// also a C-centred orthorhombic one (a_C = a+b, b_C = -a+b, c_C = c), and asking for that class has to
|
|
// return that setting even though the plain search rightly prefers the tetragonal one.
|
|
TEST_CASE("LatticeSearchForClass - tetragonal P also has a C-centred orthorhombic setting") {
|
|
const double a = 50.0, c = 120.0;
|
|
const CrystalLattice L(a, a, c, 90, 90, 90);
|
|
|
|
const auto plain = LatticeSearch(L, 1e-6);
|
|
CHECK(plain.system == gemmi::CrystalSystem::Tetragonal);
|
|
CHECK(plain.centering == 'P');
|
|
|
|
const auto ortho = LatticeSearchForClass(L, gemmi::CrystalSystem::Orthorhombic, 'C', 1e-6);
|
|
REQUIRE(ortho.has_value());
|
|
CHECK(ortho->system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(ortho->centering == 'C');
|
|
const auto uc = ortho->conventional.GetUnitCell();
|
|
// The C cell is the face diagonal on a and b, so twice the volume and a = b = a_tet * sqrt(2).
|
|
CHECK(uc.a == Catch::Approx(a * std::sqrt(2.0)).margin(1e-4));
|
|
CHECK(uc.b == Catch::Approx(a * std::sqrt(2.0)).margin(1e-4));
|
|
CHECK(uc.c == Catch::Approx(c).margin(1e-4));
|
|
CHECK(uc.alpha == Catch::Approx(90).margin(1e-4));
|
|
CHECK(uc.beta == Catch::Approx(90).margin(1e-4));
|
|
CHECK(uc.gamma == Catch::Approx(90).margin(1e-4));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearchForClass - a class the metric cannot carry is refused") {
|
|
// A general triclinic metric has no monoclinic-C setting, and an F-centred cubic lattice has no
|
|
// hexagonal-P one (its hexagonal description is R-centred).
|
|
const CrystalLattice tri(41.0, 47.0, 53.0, 71.0, 83.0, 97.0);
|
|
CHECK_FALSE(LatticeSearchForClass(tri, gemmi::CrystalSystem::Monoclinic, 'C').has_value());
|
|
|
|
const double a = 60.0;
|
|
const auto cubic_f = CrystalLattice(a, a, a, 90, 90, 90).ToPrimitive('F');
|
|
CHECK(LatticeSearch(cubic_f, 1e-6).centering == 'F');
|
|
CHECK_FALSE(LatticeSearchForClass(cubic_f, gemmi::CrystalSystem::Hexagonal, 'P').has_value());
|
|
// ... but its rhombohedral setting is there, which is what makes the refusal above a real answer
|
|
// rather than an artefact of the filter.
|
|
const auto rhomb = LatticeSearchForClass(cubic_f, gemmi::CrystalSystem::Trigonal, 'R');
|
|
REQUIRE(rhomb.has_value());
|
|
CHECK(rhomb->centering == 'R');
|
|
}
|
|
|
|
TEST_CASE("LatticeSearchForClass - asking for what the plain search found returns the same setting") {
|
|
const double a = 40.0;
|
|
const auto L = CrystalLattice(a, a, a, 90, 90, 90).ToPrimitive('I');
|
|
const auto plain = LatticeSearch(L, 1e-6);
|
|
const auto filtered = LatticeSearchForClass(L, plain.system, plain.centering, 1e-6);
|
|
REQUIRE(filtered.has_value());
|
|
CHECK(filtered->niggli_class == plain.niggli_class);
|
|
check_uc(filtered->conventional.GetUnitCell(), a, a, a, 90, 90, 90, 1e-4, 1e-4);
|
|
}
|
|
|
|
// The reduction epsilon. An exactly body-centred tetragonal lattice with c > a*sqrt(2) reduces to a
|
|
// character whose gamma is 90 EXACTLY, so the scalar product that decides the Niggli type is
|
|
// structurally zero and what a float lattice carries there is rounding. Axis-aligned that rounding
|
|
// happens to vanish - which is why the two tetragonal-I cases above pass - but every lattice the
|
|
// pipeline classifies is a refined, ROTATED one, and rotating this one about its own 4-fold is
|
|
// enough to lose the 4-fold on 38 of 60 rotations.
|
|
TEST_CASE("LatticeSearch - a body-centred tetragonal lattice keeps its 4-fold once it is rotated") {
|
|
CrystalLattice L(Coord(40, 0, 0), Coord(0, 40, 0), Coord(0, 0, 90));
|
|
L = L.ToPrimitive('I').Multiply(RotMatrix(0.3f, Coord(0, 0, 1)));
|
|
const auto res = LatticeSearch(L, 1e-6);
|
|
CHECK(res.system == gemmi::CrystalSystem::Tetragonal);
|
|
CHECK(res.centering == 'I');
|
|
}
|
|
|
|
// ITA character 43, the mI form. An ordinary centred-monoclinic crystal that happens to reduce into
|
|
// the form the table names mI - the same Bravais lattice in another setting, there is no fifteenth
|
|
// type. With that row absent the walk reaches character 44 and the centring is lost outright. The
|
|
// three monoclinic-C cases above reduce to characters 14, 39 and 14, so none of them samples it.
|
|
TEST_CASE("LatticeSearch - a centred monoclinic lattice that reduces to the mI form keeps its centring") {
|
|
const CrystalLattice L = CrystalLattice(35, 60, 30, 90, 120, 90).ToPrimitive('C');
|
|
const auto res = LatticeSearch(L, 1e-6);
|
|
CHECK(res.system == gemmi::CrystalSystem::Monoclinic);
|
|
CHECK(res.centering == 'I');
|
|
const auto uc = res.conventional.GetUnitCell();
|
|
CHECK(std::fabs(uc.alpha - 90.0) < 1e-3);
|
|
CHECK(std::fabs(uc.gamma - 90.0) < 1e-3);
|
|
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));
|
|
}
|
|
|
|
// SymmetrizeMetric is what a run falls back on when it has adopted a group on a cell no fit under
|
|
// that group produced. It has to do two things: leave the group's metric exactly satisfied, and move
|
|
// the cell as little as that requires - which means the orientation it arrives in is kept.
|
|
TEST_CASE("SymmetrizeMetric puts a cell onto the metric its group fixes") {
|
|
const auto &c2 = *gemmi::find_spacegroup_by_name("C 1 2 1");
|
|
|
|
// A C-centred monoclinic setting whose alpha is 1.5 deg off, as a freely refined metric promoted
|
|
// after integration arrives: the group cannot describe it.
|
|
const CrystalLattice off(160.0, 140.0, 90.0, 88.5, 119.9, 90.1);
|
|
const auto fixed = SymmetrizeMetric(off, c2).GetUnitCell();
|
|
CHECK(fixed.alpha == Catch::Approx(90.0).margin(1e-3));
|
|
CHECK(fixed.gamma == Catch::Approx(90.0).margin(1e-3));
|
|
// The free angle and the lengths stay where they were, to well inside the move it had to make.
|
|
CHECK(fixed.beta == Catch::Approx(119.9).margin(0.2));
|
|
CHECK(fixed.a == Catch::Approx(160.0).epsilon(2e-3));
|
|
CHECK(fixed.b == Catch::Approx(140.0).epsilon(2e-3));
|
|
CHECK(fixed.c == Catch::Approx(90.0).epsilon(2e-3));
|
|
|
|
// A cell the group already describes is not moved at all, whatever orientation it is in.
|
|
const CrystalLattice ok(80.0, 50.0, 60.0, 90.0, 105.0, 90.0);
|
|
const auto same = SymmetrizeMetric(ok, c2).GetUnitCell();
|
|
check_uc(same, 80.0, 50.0, 60.0, 90.0, 105.0, 90.0, 1e-3, 1e-3);
|
|
|
|
// ...including one that is not axis-aligned: the answer is a property of the lattice, not of the
|
|
// frame it is written in, and the orientation it came in is the orientation it goes out in.
|
|
const RotMatrix rot(0.7, Coord(1.0f, 2.0f, 3.0f).Normalize());
|
|
const CrystalLattice turned = off.Multiply(rot);
|
|
// A LEFT-handed basis keeps its hand and its angles: the three-Coord constructor would flip an
|
|
// axis and hand back the supplement of beta, which is a different cell from the one the
|
|
// reflections are indexed on - and a metric projection may not change the cell that much.
|
|
const CrystalLattice flipped = off.Multiply(gemmi::Mat33(-1, 0, 0, 0, 1, 0, 0, 0, 1));
|
|
REQUIRE(flipped.CalcVolume() < 0);
|
|
const auto left = SymmetrizeMetric(flipped, c2);
|
|
CHECK(left.CalcVolume() < 0);
|
|
CHECK(left.GetUnitCell().beta == Catch::Approx(flipped.GetUnitCell().beta).margin(1e-3));
|
|
CHECK(left.GetUnitCell().alpha == Catch::Approx(90.0).margin(1e-3));
|
|
|
|
const auto turned_fixed = SymmetrizeMetric(turned, c2);
|
|
check_uc(turned_fixed.GetUnitCell(), fixed.a, fixed.b, fixed.c, 90.0, fixed.beta, 90.0, 1e-3, 1e-3);
|
|
CHECK(turned_fixed.Vec0() * turned.Vec0()
|
|
== Catch::Approx(turned_fixed.Vec0().Length() * turned.Vec0().Length()).epsilon(1e-4));
|
|
}
|
|
|
|
TEST_CASE("LatticeSearch - a pseudo-tetragonal cell is tetragonal at the walk's tolerance and "
|
|
"orthorhombic with no length equality granted") {
|
|
// a and b 2.4 % apart: inside LATTICE_SEARCH_DIST_TOLERANCE, far outside what spot positions
|
|
// resolve. The two hypotheses rugnux carries forward are exactly these two answers.
|
|
const double a = 5.795, b = 5.933, c = 12.362;
|
|
CrystalLattice L(Coord(a, 0, 0), Coord(0, b, 0), Coord(0, 0, c));
|
|
|
|
const auto promoted = LatticeSearch(L);
|
|
CHECK(promoted.system == gemmi::CrystalSystem::Tetragonal);
|
|
CHECK(promoted.centering == 'P');
|
|
CHECK(ClassImposesLengthEquality(promoted.system));
|
|
|
|
const auto free_class = LatticeSearch(L, /*dist_tolerance=*/0.0);
|
|
CHECK(free_class.system == gemmi::CrystalSystem::Orthorhombic);
|
|
CHECK(free_class.centering == 'P');
|
|
const auto uc = free_class.conventional.GetUnitCell();
|
|
check_uc(uc, a, b, c, 90.0, 90.0, 90.0, 1e-2, 1e-2);
|
|
}
|
|
|
|
TEST_CASE("ClassImposesLengthEquality names the classes whose metric asserts a = b") {
|
|
CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Tetragonal));
|
|
CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Trigonal));
|
|
CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Hexagonal));
|
|
CHECK(ClassImposesLengthEquality(gemmi::CrystalSystem::Cubic));
|
|
CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Orthorhombic));
|
|
CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Monoclinic));
|
|
CHECK_FALSE(ClassImposesLengthEquality(gemmi::CrystalSystem::Triclinic));
|
|
}
|
|
|
|
TEST_CASE("LengthEqualityDeparture reads a = b in the class's own basis, not a stale reindex") {
|
|
// A hexagonal lattice whose conventional cell has been re-expressed in hexagonal axes while
|
|
// `reindex` still describes the C-centred orthorhombic setting the walk matched - what
|
|
// RotationIndexer does. The departure must be read through the primitive/conventional pair, not
|
|
// that matrix: the C-ortho setting has b = a*sqrt(3) by construction, so reading a against b
|
|
// there reports 53.6 % on every hexagonal lattice, whatever the crystal.
|
|
const double a = 60.0, c = 95.0;
|
|
// Primitive hexagonal cell: two equal axes at 120 deg, third perpendicular.
|
|
const CrystalLattice prim(Coord(a, 0, 0), Coord(-a / 2, a * std::sqrt(3.0) / 2, 0), Coord(0, 0, c));
|
|
|
|
LatticeSearchResult sr;
|
|
sr.system = gemmi::CrystalSystem::Hexagonal;
|
|
sr.centering = 'P';
|
|
sr.primitive_reduced = prim;
|
|
sr.conventional = prim; // already the hexagonal setting
|
|
// The stale matrix: the ortho-hexagonal C-centred setting (a, a + 2b, c).
|
|
sr.reindex = gemmi::Mat33(1, 0, 0, 1, 2, 0, 0, 0, 1);
|
|
|
|
// The free refinement of the same lattice: a and b 0.5 % apart, no more.
|
|
const double b_free = a * 1.005;
|
|
const CrystalLattice free_prim(Coord(a, 0, 0),
|
|
Coord(-b_free / 2, b_free * std::sqrt(3.0) / 2, 0),
|
|
Coord(0, 0, c));
|
|
|
|
CHECK(LengthEqualityDeparture(sr, free_prim) == Catch::Approx(0.005).margin(5e-4));
|
|
|
|
// Reading it through the stale reindex is the defect this pins: the sqrt(3) signature.
|
|
const UnitCell wrong = free_prim.Multiply(sr.reindex).GetUnitCell();
|
|
CHECK(std::fabs(wrong.a - wrong.b) / (0.5 * (wrong.a + wrong.b)) > 0.4);
|
|
}
|
|
|
|
TEST_CASE("LengthEqualityDeparture is zero for a class that asserts no length equality") {
|
|
LatticeSearchResult sr;
|
|
sr.system = gemmi::CrystalSystem::Orthorhombic;
|
|
sr.centering = 'P';
|
|
sr.primitive_reduced = CrystalLattice(Coord(10, 0, 0), Coord(0, 20, 0), Coord(0, 0, 30));
|
|
sr.conventional = sr.primitive_reduced;
|
|
CHECK(LengthEqualityDeparture(sr, sr.primitive_reduced) == 0.0);
|
|
}
|