mx/Jungfraujoch
3
2
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 81 Wiki Activity

rugnux: the model maps come from FFTW too, so the model path uses one FFT library

Browse Source 
The map-coefficient to real-space transform was still gemmi's bundled pocketfft while the
map to structure-factor direction had moved to FFTW. MapFromFPhi does that inverse with
FFTW - the same half-l XYZ layout, the same 1/V scale, NaN coefficients read as zero, the
same conjugation convention - and the output maps are now built with it. Plans are
FFTW_ESTIMATE, cached per grid size and direction, and made under the shared planner lock.

No pocketfft code is instantiated in rugnux any more (75 symbols -> 0); gemmi's fourier.hpp
is still included for its non-FFT helpers (get_size_for_hkl, get_f_phi_on_grid), so the
vendored header and its notice stay.

Tested element-wise against gemmi on odd and even grids, with arbitrary phases, negative
indices and symmetry/Friedel expansion in three space groups, plus a round trip; both new
tests fail if the conjugation is dropped. On the eight --model audit sets the merged MTZ is
byte-identical, the map coefficients are identical, the CCP4 maps agree to 3e-7 relative
(map CC 1.00000000) and every reported key is unchanged.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT
This commit is contained in:
leonarski_f
2026-09-20 18:45:17 +02:00
co-authored by Claude Opus 5
parent 8bc4c6e8a6
commit 1b414909cb
4 changed files with 131 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/ModelValidationTest.cpp
+72 -1
View File
@@ -593,10 +593,12 @@ TEST_CASE("ModelValidation_CCModelFollowsTheSignalByShell", "[ModelValidation]")
// The model path's structure factors come from FFTW; they must be gemmi's own transform to float
// precision, in the same layout, so prepare_asu_data() reads the same reflections from either.
TEST_CASE("ModelValidation_MapToFPhiMatchesGemmi", "[ModelValidation]") {
// An even grid, and an odd one on every axis (FFTW and pocketfft split odd lengths differently).
const auto size = GENERATE(std::array<int, 3>{20, 24, 30}, std::array<int, 3>{15, 21, 27});
gemmi::Grid<float> map;
map.unit_cell.set(40.0, 50.0, 60.0, 90.0, 95.0, 90.0);
map.spacegroup = gemmi::find_spacegroup_by_name("P 1");
map.set_size(20, 24, 30);
map.set_size(size[0], size[1], size[2]);
std::mt19937 rng(7);
std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
for (auto &x : map.data)
@@ -623,6 +625,75 @@ TEST_CASE("ModelValidation_MapToFPhiMatchesGemmi", "[ModelValidation]") {
CHECK(a.v[i].hkl == b.v[i].hkl);
}
// The output maps come from FFTW too; each must be gemmi's own map to float precision, point for point
// in the same layout. Coefficients with arbitrary phases, negative indices, expanded by symmetry and
// Friedel, on grids that are odd along u and v (w is even by construction of a half-l grid).
TEST_CASE("ModelValidation_MapFromFPhiMatchesGemmi", "[ModelValidation]") {
const char *sg_name = GENERATE("P 1", "P 1 21 1", "P 21 21 21");
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name(sg_name);
REQUIRE(sg != nullptr);
gemmi::AsuData<std::complex<float>> coef;
coef.unit_cell_.set(31.0, 43.0, 57.0, 90.0, sg->number == 4 ? 104.0 : 90.0, 90.0);
coef.spacegroup_ = sg;
const gemmi::ReciprocalAsu asu(sg);
const gemmi::GroupOps gops = sg->operations();
std::mt19937 rng(11);
std::uniform_real_distribution<float> amp(0.1f, 10.0f), phase(-3.14159f, 3.14159f);
for (int h = -7; h <= 7; h++)
for (int k = -9; k <= 9; k++)
for (int l = -11; l <= 11; l++) {
const gemmi::Op::Miller hkl{{h, k, l}};
if ((h == 0 && k == 0 && l == 0) || !asu.is_in(hkl) || gops.is_systematically_absent(hkl))
continue;
coef.v.push_back({hkl, std::polar(amp(rng), phase(rng))});
}
REQUIRE(coef.v.size() > 500);
// P 1 takes an odd grid on u and v; the screw axes need even factors.
const std::array<int, 3> size = sg->number == 1 ? std::array<int, 3>{17, 21, 26}
: std::array<int, 3>{18, 24, 26};
gemmi::FPhiGrid<float> grid = gemmi::get_f_phi_on_grid<float>(coef, size, true);
grid.data[grid.index_n(2, -3, 4)] = std::complex<float>(1.0f, NAN); // a missing coefficient
const gemmi::Grid<float> ours = MapFromFPhi(grid);
const gemmi::Grid<float> ref = gemmi::transform_f_phi_grid_to_map(gemmi::FPhiGrid<float>(grid));
REQUIRE(ours.nu == ref.nu);
REQUIRE(ours.nv == ref.nv);
REQUIRE(ours.nw == ref.nw);
REQUIRE(ours.axis_order == ref.axis_order);
REQUIRE(ours.spacegroup == ref.spacegroup);
REQUIRE(ours.data.size() == ref.data.size());
double largest = 0, worst = 0;
for (size_t i = 0; i < ref.data.size(); i++) {
REQUIRE(std::isfinite(ours.data[i]));
largest = std::max(largest, static_cast<double>(std::abs(ref.data[i])));
worst = std::max(worst, static_cast<double>(std::abs(ours.data[i] - ref.data[i])));
}
CHECK(largest > 0);
CHECK(worst <= 1e-5 * largest);
}
// Map -> coefficients -> map is the identity (the V/N and 1/V scales cancel the unnormalised
// transforms), on a grid odd along u and v.
TEST_CASE("ModelValidation_ModelFFTRoundTrip", "[ModelValidation]") {
gemmi::Grid<float> map;
map.unit_cell.set(35.0, 45.0, 55.0, 80.0, 95.0, 105.0);
map.spacegroup = gemmi::find_spacegroup_by_name("P 1");
map.set_size(15, 21, 28);
std::mt19937 rng(3);
std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
for (auto &x : map.data)
x = dist(rng);
const gemmi::Grid<float> back = MapFromFPhi(MapToFPhi(map));
REQUIRE(back.nu == map.nu);
REQUIRE(back.nv == map.nv);
REQUIRE(back.nw == map.nw);
double worst = 0;
for (size_t i = 0; i < map.data.size(); i++)
worst = std::max(worst, static_cast<double>(std::abs(back.data[i] - map.data[i])));
CHECK(worst <= 1e-5);
}
// The Jacobian's six columns are evaluated in parallel; the placement must be the serial one, bit for
// bit.
TEST_CASE("ModelValidation_RigidBodySameOnAnyNumberOfThreads", "[ModelValidation]") {