mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions 2 Packages Projects Releases 73 Wiki Activity
Files
6491345feaa8290ee30f30217903ad110bd4096b
Jungfraujoch/tests/CrystalLatticeTest.cpp
T
leonarski_fandClaude Opus 5 4e8db41785 indexing: refuse a coplanar candidate, and search the plane normal when the shortlist is flat 
Three related changes to the FFT candidate path, batteried together because they touch the same
function.

A COPLANAR CANDIDATE REACHED REFINEMENT. ReduceResults filtered triples on lengths and angles only -
the 30-150 degree bound admits any flat combination - and there was no volume test. On one dataset 41
of 5535 candidates had |V|/abc below 0.05, with a clean decade gap to the next, and three of them
reached the optimizer. UnitCell is float, and for a cell that flat the metric determinant is around
1.5e-7, so float32 gets its sign wrong 19% of the time where float64 never does. The guard against a
negative argument to sqrt then CREATES the singularity it was meant to prevent: it puts c in the a-b
plane, the reciprocal volume is 1/0, and the residual is 0 times infinity. Ceres reported a
not-a-number Jacobian and wrote several hundred lines of solver output per failed solve.

VolumeFraction() is |V|/(|a||b||c|), rejected below 0.02 - about 1.1 degrees off flat, ten times below
the flattest real candidate observed and a thousand times above where float loses the sign. It is
enforced at the producer and at the two optimizer entry points. Note the existing sanity checks use
ABSOLUTE volume, which a 320 cubic-angstrom flat cell passes. The same reciprocal-volume division is
now guarded at the two remaining sites that share the pattern.

A SHORTLIST CONFINED TO ONE PLANE cannot close a cell, and the row it is missing is the plane normal.
That is detected from the scatter-matrix eigenvalue ratio - measured, degenerate clouds score 2e-5 to
3.3e-4 against 0.026 or more for every non-degenerate one, a factor of eighty - and one further
transform is spent with the same direction count inside a three-degree cap about the normal, so the
plan and buffers are untouched. More directions cannot substitute: at the exact true direction the
long axis ranks 1422 of 16384 by prominence while the shortlist cut is four times higher. Ranking, not
sampling, is the obstacle. A four-fold denser grid was measured and rejected - it reaches the same
answer to three decimal places and takes a run from 2.5 to 8 GB of device memory.

fft_min_unit_cell_A is reachable as --fft-min-unit-cell and is lowered automatically by -C, mirroring
how the maximum is already raised. The default of 10 is unchanged: a lower floor admits spurious
sub-cells on protein data, and over 73 protein runs the floor was never lowered while the sibling
maximum did fire twice, so the path is live and correctly inert.

Corpus of 93 datasets, both arms, one build: 72 bit-identical on report content and p.hkl checksum, 13
failing identically, and the count of working datasets rises by one. The volume guard fires on 58 of
93 and 47 of those stay bit-identical - it fires constantly and almost never changes an answer, which
is what it should do. Solver chatter falls from 919 lines across three datasets to none. The cap
fires on 4 of 93, none of them in the in-house or private arms.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01Lc5JG6kJqZoCWaoZ43JGTW
2026-08-29 20:31:20 +02:00

210 lines
8.1 KiB
C++
Raw Blame History

// SPDX-FileCopyrightText: 2024 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"
TEST_CASE("CrystalLattice") {
CrystalLattice l(50,60,80, 90, 90, 90);
REQUIRE(l.Vec0().Length() == Catch::Approx(50));
REQUIRE(l.Vec1().Length() == Catch::Approx(60));
REQUIRE(l.Vec2().Length() == Catch::Approx(80));
REQUIRE(angle_deg(l.Vec0(), l.Vec2()) == 90);
REQUIRE(angle_deg(l.Vec0(), l.Vec1()) == 90);
REQUIRE(angle_deg(l.Vec1(), l.Vec2()) == 90);
auto uc0 = l.GetUnitCell();
REQUIRE(uc0.a == Catch::Approx(50));
REQUIRE(uc0.b == Catch::Approx(60));
REQUIRE(uc0.c == Catch::Approx(80));
REQUIRE(uc0.alpha == Catch::Approx(90));
REQUIRE(uc0.beta == Catch::Approx(90));
REQUIRE(uc0.gamma == Catch::Approx(90));
l = CrystalLattice(30, 40, 70, 90, 95, 90);
REQUIRE(l.Vec0().Length() == Catch::Approx(30));
REQUIRE(l.Vec1().Length() == Catch::Approx(40));
REQUIRE(l.Vec2().Length() == Catch::Approx(70));
REQUIRE(angle_deg(l.Vec0(), l.Vec2()) == 95);
REQUIRE(angle_deg(l.Vec0(), l.Vec1()) == 90);
REQUIRE(angle_deg(l.Vec1(), l.Vec2()) == 90);
auto uc1 = l.GetUnitCell();
REQUIRE(uc1.a == Catch::Approx(30));
REQUIRE(uc1.b == Catch::Approx(40));
REQUIRE(uc1.c == Catch::Approx(70));
REQUIRE(uc1.alpha == Catch::Approx(90));
REQUIRE(uc1.beta == Catch::Approx(95));
REQUIRE(uc1.gamma == Catch::Approx(90));
l = CrystalLattice(45, 45, 70, 90, 90, 120);
REQUIRE(l.Vec0().Length() == Catch::Approx(45));
REQUIRE(l.Vec1().Length() == Catch::Approx(45));
REQUIRE(l.Vec2().Length() == Catch::Approx(70));
REQUIRE(angle_deg(l.Vec0(), l.Vec2()) == Catch::Approx(90));
REQUIRE(angle_deg(l.Vec0(), l.Vec1()) == Catch::Approx(120));
REQUIRE(angle_deg(l.Vec1(), l.Vec2()) == Catch::Approx(90));
auto uc2 = l.GetUnitCell();
REQUIRE(uc2.a == Catch::Approx(45));
REQUIRE(uc2.b == Catch::Approx(45));
REQUIRE(uc2.c == Catch::Approx(70));
REQUIRE(uc2.alpha == Catch::Approx(90));
REQUIRE(uc2.beta == Catch::Approx(90));
REQUIRE(uc2.gamma == Catch::Approx(120));
}
TEST_CASE("CrystalLattice_Sort") {
CrystalLattice l(80,60,50, 120, 90, 90);
l.Sort();
REQUIRE(l.Vec0().Length() == Catch::Approx(50));
REQUIRE(l.Vec1().Length() == Catch::Approx(60));
REQUIRE(l.Vec2().Length() == Catch::Approx(80));
REQUIRE(angle_deg(l.Vec0(), l.Vec2()) == Catch::Approx(90));
REQUIRE(angle_deg(l.Vec0(), l.Vec1()) == Catch::Approx(120));
REQUIRE(angle_deg(l.Vec1(), l.Vec2()) == Catch::Approx(90));
}
TEST_CASE("CrystalLattice_ReorderMonoclinic") {
std::vector<CrystalLattice> latt = {
{60, 85, 70, 90, 70, 90},
{60, 85, 70, 90, 110, 90},
};
for (const auto &l_in :latt) {
CrystalLattice l = l_in;
l.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(110.0));
}
}
TEST_CASE("CrystalLattice_ReorderTetragonal") {
CrystalLattice l(40,60,60, 90, 90, 90);
l.ReorderABEqual();
REQUIRE(l.Vec0().Length() == Catch::Approx(60));
REQUIRE(l.Vec1().Length() == Catch::Approx(60));
REQUIRE(l.Vec2().Length() == Catch::Approx(40));
}
TEST_CASE("CrystalLattice_ReorderTetragonal_ForHexagonal") {
CrystalLattice l(40,60,60, 120, 90, 90);
l.ReorderABEqual();
REQUIRE(l.Vec0().Length() == Catch::Approx(60));
REQUIRE(l.Vec1().Length() == Catch::Approx(60));
REQUIRE(l.Vec2().Length() == Catch::Approx(40));
}
TEST_CASE("CrystalLattice_Handedness") {
CrystalLattice l(Coord(1,0,0), Coord(0,1,0), Coord(0,0,-1));
REQUIRE(l.Vec0().x == Catch::Approx(1));
REQUIRE(l.Vec1().y == Catch::Approx(1));
REQUIRE(l.Vec2().z == Catch::Approx(1));
}
TEST_CASE("CrystalLattice_Volume") {
CrystalLattice l(50, 60, 80, 90, 90, 90);
REQUIRE(l.CalcVolume() == 50 * 60 * 80);
CrystalLattice l2(50,60,80, 90, 120, 90);
float sin120 = std::sqrt(3) / 2;
REQUIRE(l2.CalcVolume() == Catch::Approx(50 * 60 * 80 * sin120));
}
TEST_CASE("CrystalLattice_VolumeFraction") {
CrystalLattice ortho(50, 60, 80, 90, 90, 90);
REQUIRE(ortho.VolumeFraction() == Catch::Approx(1.0));
// All three angles at 60 deg is about as oblique as a reduced cell gets, and is still an order
// of magnitude clear of MIN_BASIS_VOLUME_FRACTION.
CrystalLattice oblique(50, 50, 50, 60, 60, 60);
REQUIRE(oblique.VolumeFraction() == Catch::Approx(std::sqrt(0.5)).margin(1e-4));
REQUIRE(oblique.VolumeFraction() > MIN_BASIS_VOLUME_FRACTION);
// Three coplanar rows: ordinary lengths, ordinary angles, no volume. This is what an indexing
// candidate built from directions that all lie in one dense reciprocal plane looks like, and
// what nothing downstream can refine - 1/V is not finite.
CrystalLattice flat(Coord(50, 0, 0), Coord(0, 60, 0), Coord(30, 40, 0));
REQUIRE(flat.VolumeFraction() < MIN_BASIS_VOLUME_FRACTION);
CrystalLattice degenerate(Coord(0, 0, 0), Coord(0, 60, 0), Coord(0, 0, 80));
REQUIRE(degenerate.VolumeFraction() == 0.0f);
}
TEST_CASE("CrystalLattice_Recip") {
CrystalLattice l(50,60,80, 90, 90, 90);
REQUIRE(l.Astar().Length() == Catch::Approx(1/50.0));
REQUIRE(l.Astar().x == Catch::Approx(1/50.0));
REQUIRE(l.Bstar().Length() == Catch::Approx(1/60.0));
REQUIRE(l.Bstar().y == Catch::Approx(1/60.0));
REQUIRE(l.Cstar().Length() == Catch::Approx(1/80.0));
REQUIRE(l.Cstar().z == Catch::Approx(1/80.0));
}
TEST_CASE("CrystalLattice_ToPrimitive") {
// Conventional cubic (take as I-conventional basis): a=b=c=50, orthogonal
CrystalLattice conv(50, 50, 50, 90, 90, 90);
const float Vconv = conv.CalcVolume();
REQUIRE(Vconv == Catch::Approx(50.f * 50.f * 50.f).margin(1e-4f));
CrystalLattice prim = conv.ToPrimitive('I');
// Volume halves (index 2)
CHECK(prim.CalcVolume() == Catch::Approx(Vconv * 0.5f).margin(1e-4f));
// bcc primitive lengths and angles
auto uc_p = prim.GetUnitCell();
const float a0 = 50.f;
const float expected_len = a0 * std::sqrt(3.f) * 0.5f; // a*sqrt(3)/2
const float expected_angle = 109.4712206f; // arccos(-1/3) in degrees
CHECK(uc_p.a == Catch::Approx(expected_len).margin(1e-3f));
CHECK(uc_p.b == Catch::Approx(expected_len).margin(1e-3f));
CHECK(uc_p.c == Catch::Approx(expected_len).margin(1e-3f));
CHECK(uc_p.alpha == Catch::Approx(expected_angle).margin(1e-3f));
CHECK(uc_p.beta == Catch::Approx(expected_angle).margin(1e-3f));
CHECK(uc_p.gamma == Catch::Approx(expected_angle).margin(1e-3f));
}
TEST_CASE("CrystalLattice_FromPrimitive") {
// Conventional cubic (take as I-conventional basis): a=b=c=50, orthogonal
const float a0 = 50.f;
const float expected_len = a0 * std::sqrt(3.f) * 0.5f; // a*sqrt(3)/2
const float expected_angle = 109.4712206f; // arccos(-1/3) in degrees
CrystalLattice conv(expected_len, expected_len, expected_len,
expected_angle, expected_angle, expected_angle);
auto prim = conv.FromPrimitive('I');
auto uc_c = prim.GetUnitCell();
// Back to conventional cube with original volume
CHECK(prim.CalcVolume() == Catch::Approx(a0 * a0 * a0).margin(1e-4f));
CHECK(uc_c.a == Catch::Approx(a0).margin(1e-3f));
CHECK(uc_c.b == Catch::Approx(a0).margin(1e-3f));
CHECK(uc_c.c == Catch::Approx(a0).margin(1e-3f));
CHECK(uc_c.alpha == Catch::Approx(90.f).margin(1e-3f));
CHECK(uc_c.beta == Catch::Approx(90.f).margin(1e-3f));
CHECK(uc_c.gamma == Catch::Approx(90.f).margin(1e-3f));
}
TEST_CASE("CrystalLattice_NiggliReduce") {
Coord a = {3,0,0};
Coord b = {0,4,0};
Coord c = {0,0,5};
CrystalLattice l(a + 2 * b, b, c - 5 * b);
auto red_uc = l.NiggliReduce().GetUnitCell();
CHECK(red_uc.a == Catch::Approx(3.0));
CHECK(red_uc.b == Catch::Approx(4.0));
CHECK(red_uc.c == Catch::Approx(5.0));
CHECK(red_uc.alpha == Catch::Approx(90.0));
CHECK(red_uc.beta == Catch::Approx(90.0));
CHECK(red_uc.gamma == Catch::Approx(90.0));
}
Reference in New Issue View Git Blame Copy Permalink