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

rugnux --model: grid the rigid body's density and solvent mask on all threads

Browse Source 
The rigid-body target re-grids the model at every evaluation (gemmi's
put_model_density_on_grid with the refmac blur, then the Refmac solvent
mask), and that gridding was ~80% of an evaluation on low-symmetry cells.
The central evaluation ran it on one thread and the six Jacobian columns
on one pool worker each, so most of a 32-thread machine sat idle through
the real fit.

New rugnux/ModelGrid.{h,cpp} reimplements the two gemmi routines so that
they are bit-identical to gemmi on any number of threads:

- Atoms are put on the grid per w plane: each plane is one task and walks
  every atom whose box reaches it, in model order, visiting only its own
  points (a copy of gemmi's do_use_points_in_box restricted to one plane).
  Every grid point therefore receives the same additions in the same
  order as in gemmi's serial loop. Per-atom coefficients and radius are
  computed once up front, and each atom's density function is copied to a
  local so the compiler keeps it in registers.
- Symmetrization runs over orbits: gemmi reduces each orbit into its
  lowest index, so the leaders are the points with no lower mate. They
  depend only on the grid size and group, so RefineRigidBody finds them
  once per zone; each orbit is then reduced by one thread with gemmi's
  operands in gemmi's order. Orbits share no points.
- The solvent mask uses the same two passes (setting points to 0 is order
  independent anyway); gemmi's own island removal and shrink follow.
  vendored gemmi is untouched.

The six Jacobian columns now run on std::async threads, not on the pool,
so each column's gridding can spread over the pool (a pass reached from a
pool worker runs inline). With one thread they run deferred, serially.

Evidence: the grids are memcmp-identical to gemmi's at nt 1/5/32 on 41
deposited models from the battery's PDB cache plus the 8 battery sets
(23 of them with anisotropic atoms; P1 up to F4132, R3/R32, I and C
centring) at 6/4.5/3.5 A; new test ModelValidation_ParallelGriddingMatchesGemmi
covers P1, C2, P212121, I23 and F4132 with iso and aniso atoms. A
standalone RefineRigidBody benchmark (synthetic |F| from the model,
displaced model) gives identical evaluations/angle/shift/coordinates to
the unchanged code: 5lzl 24-33 s -> 12 s real fit on a loaded machine;
null (9 concurrent replicates) unchanged within noise.

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C
This commit is contained in:
leonarski_f
2026-09-28 02:02:28 +02:00
co-authored by Claude Opus 5.5
parent 530c4feb90
commit 158fb0edbf
6 changed files with 392 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/ModelValidationTest.cpp
+69
View File
@@ -5,6 +5,7 @@
#include <cmath>
#include <cstdio>
#include <cstring>
#include <filesystem>
#include <fstream>
#include <map>
@@ -13,9 +14,11 @@
#include <gemmi/mmread_gz.hpp>
#include <gemmi/fourier.hpp>
#include <gemmi/solmask.hpp>
#include "../common/Logger.h"
#include "../rugnux/ModelFFT.h"
#include "../rugnux/ModelGrid.h"
#include "../rugnux/ModelValidation.h"
#include "../rugnux/RigidBodyRefine.h"
#include "../rugnux/SigmaA.h"
@@ -811,3 +814,69 @@ TEST_CASE("ModelValidation_RigidBodySameOnAnyNumberOfThreads", "[ModelValidation
std::filesystem::remove(path);
}
// The rigid body puts each probe placement on the grid with its own parallel copy of gemmi's gridding
// (density, solvent mask and their symmetrization). It must give gemmi's grids bit for bit, on any
// number of threads, for isotropic and anisotropic atoms and for groups with and without centring.
TEST_CASE("ModelValidation_ParallelGriddingMatchesGemmi", "[ModelValidation]") {
const char *crysts[] = {
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 1 1\n",
"CRYST1 40.000 50.000 60.000 90.00 100.00 90.00 C 1 2 1 4\n",
"CRYST1 40.000 50.000 60.000 90.00 90.00 90.00 P 21 21 21 4\n",
"CRYST1 60.000 60.000 60.000 90.00 90.00 90.00 I 2 3 24\n",
"CRYST1 80.000 80.000 80.000 90.00 90.00 90.00 F 41 3 2 96\n",
};
for (const char *cryst : crysts) {
const auto path = WriteTemp("parallel_gridding_test.pdb", ClusterPdb(cryst).c_str());
gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
std::filesystem::remove(path);
const gemmi::SpaceGroup *sg = st.find_spacegroup();
REQUIRE(sg != nullptr);
int i = 0;
for (gemmi::Chain &ch : st.models[0].chains)
for (gemmi::Residue &r : ch.residues)
for (gemmi::Atom &a : r.atoms)
if (i++ % 3 == 0)
a.aniso = {0.30f, 0.25f, 0.20f, 0.02f, -0.01f, 0.03f};
const auto same = [](const gemmi::Grid<float> &a, const gemmi::Grid<float> &b) {
return a.data.size() == b.data.size() &&
std::memcmp(a.data.data(), b.data.data(), a.data.size() * sizeof(float)) == 0;
};
for (double d_min : {6.0, 3.5}) {
gemmi::DensityCalculator<gemmi::IT92<float>, float> gemmi_dc;
gemmi_dc.d_min = d_min;
gemmi_dc.rate = 1.5;
gemmi_dc.grid.unit_cell = st.cell;
gemmi_dc.grid.spacegroup = sg;
gemmi_dc.set_refmac_compatible_blur(st.models[0]);
gemmi_dc.put_model_density_on_grid(st.models[0]);
gemmi::Grid<float> gemmi_mask;
gemmi_mask.unit_cell = st.cell;
gemmi_mask.spacegroup = sg;
gemmi_mask.set_size_from_spacing(gemmi_dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
gemmi::SolventMasker(gemmi::AtomicRadiiSet::Refmac).put_mask_on_grid(gemmi_mask, st.models[0]);
for (size_t nthreads : {1, 4}) {
const std::vector<size_t> leaders = OrbitLeaders(gemmi_mask, nthreads);
gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
dc.d_min = d_min;
dc.rate = 1.5;
dc.grid.unit_cell = st.cell;
dc.grid.spacegroup = sg;
dc.set_refmac_compatible_blur(st.models[0]);
PutModelDensityOnGrid(dc, st.models[0], leaders, nthreads);
CHECK(same(dc.grid, gemmi_dc.grid));
gemmi::Grid<float> mask;
mask.unit_cell = st.cell;
mask.spacegroup = sg;
mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
PutMaskOnGrid(mask, st.models[0], leaders, nthreads);
CHECK(same(mask, gemmi_mask));
}
}
}
}