From 158fb0edbf12a32298319a04fecc55d77e25208a Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 28 Sep 2026 02:02:28 +0200 Subject: [PATCH] rugnux --model: grid the rigid body's density and solvent mask on all threads 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) Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C --- docs/CHANGELOG.md | 1 + rugnux/CMakeLists.txt | 2 + rugnux/ModelGrid.cpp | 245 ++++++++++++++++++++++++++++++++++ rugnux/ModelGrid.h | 29 ++++ rugnux/RigidBodyRefine.cpp | 68 +++++++--- tests/ModelValidationTest.cpp | 69 ++++++++++ 6 files changed, 392 insertions(+), 22 deletions(-) create mode 100644 rugnux/ModelGrid.cpp create mode 100644 rugnux/ModelGrid.h diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index b9ab621ee..da57c124d 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -8,6 +8,7 @@ * jfjoch_viewer keeps a separate preferred dataset-info plot for grid scans, where "Spots + background" means the spot count. * jfjoch_viewer's dark theme no longer leaves navy buttons, red warnings and chart guide lines at their light-theme colours. * Rugnux adds beam-stop holder arms that let part of the beam through to the beam-stop mask, without changing the mask of a sweep that has none. +* Rugnux's `--model` rigid-body refinement computes the model's density and solvent mask on all threads, with results unchanged. ### 1.0.0-rc.173 diff --git a/rugnux/CMakeLists.txt b/rugnux/CMakeLists.txt index b9a9c4835..ad785e24a 100644 --- a/rugnux/CMakeLists.txt +++ b/rugnux/CMakeLists.txt @@ -10,6 +10,8 @@ ADD_LIBRARY(JFJochRugnux STATIC RugnuxCommandLine.h ModelFFT.cpp ModelFFT.h + ModelGrid.cpp + ModelGrid.h ModelScaling.cpp ModelScaling.h ModelValidation.cpp diff --git a/rugnux/ModelGrid.cpp b/rugnux/ModelGrid.cpp new file mode 100644 index 000000000..125c987f2 --- /dev/null +++ b/rugnux/ModelGrid.cpp @@ -0,0 +1,245 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include "ModelGrid.h" + +#include +#include +#include +#include + +#include "gemmi/solmask.hpp" // SolventMasker + +#include "../common/ParallelFor.h" + +namespace { + +using Table = gemmi::IT92; + +// The box of grid points gemmi's Grid::use_points_in_box() walks around one atom. +struct AtomBox { + gemmi::Fractional fpos; + int du = 0, dv = 0, dw = 0; + double radius = 0; +}; + +// The box of use_points_around(fpos, radius, ..., false), sized and clamped as gemmi does it. +AtomBox MakeAtomBox(const gemmi::Grid &grid, const gemmi::Fractional &fpos, double radius) { + AtomBox box; + box.fpos = fpos; + box.radius = radius; + box.du = std::min(static_cast(std::ceil(radius / grid.spacing[0])), grid.nu - 1); + box.dv = std::min(static_cast(std::ceil(radius / grid.spacing[1])), grid.nv - 1); + box.dw = std::min(static_cast(std::ceil(radius / grid.spacing[2])), grid.nw - 1); + return box; +} + +// gemmi's use_points_in_box() run over a list of atoms (grid.hpp, do_use_points_in_box), in +// parallel over the grid's w planes: a plane is one task, and it walks every atom whose box reaches +// it, in list order, visiting only that plane's points. Every point is so handed exactly the calls +// gemmi's serial loop over the atoms hands it, in the same order, and the grid comes out bit for bit +// the same whatever the thread count. func_for(i) gives atom i's func(point, d2, delta), taken as a +// local copy so that the compiler can hold it in registers while the atom's points are written. +template +void UseAtomBoxesByPlane(gemmi::Grid &grid, const std::vector &boxes, size_t nthreads, + FuncFor func_for) { + const int nu = grid.nu, nv = grid.nv, nw = grid.nw; + std::vector> plane_atoms(nw); + for (int i = 0; i < static_cast(boxes.size()); i++) { + const int w0 = gemmi::iround(boxes[i].fpos.z * nw); + for (int w = w0 - boxes[i].dw; w <= w0 + boxes[i].dw; w++) { + std::vector &atoms = plane_atoms[gemmi::modulo(w, nw)]; + if (atoms.empty() || atoms.back() != i) + atoms.push_back(i); + } + } + ParallelFor(nw, nthreads, [&](int plane) { + for (int i : plane_atoms[plane]) { + const AtomBox &b = boxes[i]; + const auto func = func_for(i); + const double max_dist_sq = b.radius * b.radius; + const gemmi::Fractional nctr(b.fpos.x * nu, b.fpos.y * nv, b.fpos.z * nw); + const int u_lo = gemmi::iround(nctr.x) - b.du, u_hi = gemmi::iround(nctr.x) + b.du; + const int v_lo = gemmi::iround(nctr.y) - b.dv, v_hi = gemmi::iround(nctr.y) + b.dv; + const int w_lo = gemmi::iround(nctr.z) - b.dw, w_hi = gemmi::iround(nctr.z) + b.dw; + const int u_0 = gemmi::modulo(u_lo, nu), v_0 = gemmi::modulo(v_lo, nv); + gemmi::Fractional fdelta(nctr.x - u_lo, 0, 0); + // The box's w that fall on this plane, ascending as gemmi walks them: one, unless the box is + // wider than the cell. + for (int w = w_lo + gemmi::modulo(plane - w_lo, nw); w <= w_hi; w += nw) { + fdelta.z = nctr.z - w; + for (int v = v_lo, v_ = v_0; v <= v_hi; ++v, v_ = (v_ + 1 == nv ? 0 : v_ + 1)) { + fdelta.y = nctr.y - v; + gemmi::Position delta(grid.orth_n.multiply(fdelta)); + const double dist_sq0 = gemmi::sq(delta.y) + gemmi::sq(delta.z); + if (dist_sq0 > max_dist_sq) + continue; + float *t = &grid.data[grid.index_q(u_0, v_, plane)]; + for (int u = u_lo, u_ = u_0;;) { + const double dist_sq = dist_sq0 + gemmi::sq(delta.x); + if (!(dist_sq > max_dist_sq)) + func(*t, dist_sq, delta); + if (u >= u_hi) + break; + ++u; + ++u_; + ++t; + if (u_ == nu) { + u_ = 0; + t -= nu; + } + delta.x -= grid.orth_n.a11; + } + } + } + } + }); +} + +// gemmi's Grid::symmetrize() in parallel over the orbits found by OrbitLeaders(): each orbit is reduced +// by one thread, into its leader, with gemmi's operands in gemmi's order, and the value written to its +// leader and mates as gemmi writes it. Orbits share no points, so the grid is gemmi's, bit for bit. +template +void SymmetrizeOrbits(gemmi::Grid &grid, const std::vector &leaders, size_t nthreads, + Func func) { + const std::vector ops = grid.get_scaled_ops_except_id(); + const size_t nu = grid.nu, nv = grid.nv; + ParallelChunks(static_cast(leaders.size()), nthreads, [&](int lo, int hi) { + std::vector mates(ops.size()); + for (int i = lo; i < hi; i++) { + const size_t idx = leaders[i]; + const int u = static_cast(idx % nu), v = static_cast(idx / nu % nv), + w = static_cast(idx / (nu * nv)); + for (size_t k = 0; k < ops.size(); ++k) { + const std::array t = ops[k].apply(u, v, w); + mates[k] = grid.index_n(t[0], t[1], t[2]); + } + float value = grid.data[idx]; + for (size_t k : mates) + value = func(value, grid.data[k]); + grid.data[idx] = value; + for (size_t k : mates) + grid.data[k] = value; + } + }); +} + +std::vector ModelAtoms(const gemmi::Model &model) { + std::vector atoms; + for (const gemmi::Chain &ch : model.chains) + for (const gemmi::Residue &r : ch.residues) + for (const gemmi::Atom &a : r.atoms) + atoms.push_back(&a); + return atoms; +} + +// One atom's density, added to a grid point as gemmi's do_add_atom_density_to_grid() adds it. +struct AtomDensity { + using CReal = Table::Coef::coef_type; + using IsoSum = decltype(std::declval().precalculate_density_iso(CReal(), CReal())); + using AnisoSum = decltype(std::declval().precalculate_density_aniso_b( + gemmi::SMat33(), CReal())); + + bool is_aniso = false; + float occ = 0; + IsoSum iso{}; + AnisoSum aniso{}; + + void operator()(float &point, double r2, const gemmi::Position &delta) const { + if (!is_aniso) + point += float(occ * iso.calculate((CReal)r2)); + else + point += float(occ * aniso.calculate(delta)); + } +}; + +} // namespace + +// The points gemmi's Grid::symmetrize() (grid.hpp) reduces the grid's orbits into, ascending. gemmi +// walks the points by index and takes the orbit of every point it has not visited yet, so it reduces +// each orbit into the orbit's lowest index - the point none of whose mates has a lower one. They depend +// on the grid's size and space group only, so a zone finds them once for all its evaluations. +std::vector OrbitLeaders(const gemmi::Grid &grid, size_t nthreads) { + const std::vector ops = grid.get_scaled_ops_except_id(); + if (ops.empty()) + return {}; // P1: nothing to symmetrize + std::vector> plane_leaders(grid.nw); + ParallelFor(grid.nw, nthreads, [&](int w) { + for (int v = 0; v < grid.nv; ++v) + for (int u = 0; u < grid.nu; ++u) { + const size_t idx = grid.index_q(u, v, w); + bool lowest = true; + for (size_t k = 0; k < ops.size() && lowest; ++k) { + const std::array t = ops[k].apply(u, v, w); + lowest = grid.index_n(t[0], t[1], t[2]) >= idx; + } + if (lowest) + plane_leaders[w].push_back(idx); + } + }); + std::vector leaders; + for (const std::vector &l : plane_leaders) + leaders.insert(leaders.end(), l.begin(), l.end()); + return leaders; +} + +// DensityCalculator::put_model_density_on_grid() (gemmi dencalc.hpp, do_add_atom_density_to_grid) +// with the atoms and the symmetry spread over the threads as above - the same grid, bit for bit. Each +// atom's density coefficients and radius are worked out once, up front, rather than once per plane. +void PutModelDensityOnGrid(gemmi::DensityCalculator &dc, const gemmi::Model &model, + const std::vector &orbit_leaders, size_t nthreads) { + using CReal = AtomDensity::CReal; + dc.initialize_grid(); + const std::vector atoms = ModelAtoms(model); + const int n = static_cast(atoms.size()); + std::vector boxes(n); + std::vector density(n); + ParallelChunks(n, nthreads, [&](int lo, int hi) { + for (int i = lo; i < hi; i++) { + const gemmi::Atom &atom = *atoms[i]; + const auto &coef = Table::get(atom.element, atom.charge, atom.serial); + const float addend = dc.addends.get(atom.element); + const gemmi::Fractional fpos = dc.grid.unit_cell.fractionalize(atom.pos); + AtomDensity &d = density[i]; + d.occ = atom.occ; + d.is_aniso = atom.aniso.nonzero(); + if (!d.is_aniso) { + const CReal b = static_cast(atom.b_iso + dc.blur); + d.iso = coef.precalculate_density_iso(b, addend); + boxes[i] = MakeAtomBox(dc.grid, fpos, dc.estimate_radius(d.iso, b)); + } else { + const auto aniso_b = atom.aniso.scaled(CReal(gemmi::u_to_b())).added_kI(CReal(dc.blur)); + const CReal b_max = std::max(std::max(aniso_b.u11, aniso_b.u22), aniso_b.u33); + const double radius = dc.estimate_radius(coef.precalculate_density_iso(b_max, addend), b_max); + d.aniso = coef.precalculate_density_aniso_b(aniso_b, addend); + boxes[i] = MakeAtomBox(dc.grid, fpos, radius); + } + } + }); + UseAtomBoxesByPlane(dc.grid, boxes, nthreads, [&](int i) { return density[i]; }); + SymmetrizeOrbits(dc.grid, orbit_leaders, nthreads, [](float a, float b) { return a + b; }); +} + +// SolventMasker(AtomicRadiiSet::Refmac).put_mask_on_grid() (gemmi solmask.hpp) with the atoms and the +// symmetry spread over the threads as above - the same mask, bit for bit. The island removal and the +// shrink that follow are gemmi's own. +void PutMaskOnGrid(gemmi::Grid &grid, const gemmi::Model &model, const std::vector &orbit_leaders, + size_t nthreads) { + const gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac); + masker.clear(grid); + std::vector boxes; + for (const gemmi::Atom *atom : ModelAtoms(model)) { + if ((masker.ignore_hydrogen && atom->is_hydrogen()) || + (masker.ignore_zero_occupancy_atoms && atom->occ <= 0)) + continue; + double r = masker.constant_r + masker.rprobe; + r += gemmi::refmac_radius_for_bulk_solvent(atom->element.elem); + boxes.push_back(MakeAtomBox(grid, grid.unit_cell.fractionalize(atom->pos), r)); + } + UseAtomBoxesByPlane(grid, boxes, nthreads, [](int) { + return [](float &point, double, const gemmi::Position &) { point = 0.f; }; + }); + SymmetrizeOrbits(grid, orbit_leaders, nthreads, [](float a, float b) { return a < b ? a : b; }); + masker.remove_islands(grid); + masker.shrink(grid); +} diff --git a/rugnux/ModelGrid.h b/rugnux/ModelGrid.h new file mode 100644 index 000000000..4998f2226 --- /dev/null +++ b/rugnux/ModelGrid.h @@ -0,0 +1,29 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include + +#include "gemmi/dencalc.hpp" // DensityCalculator +#include "gemmi/grid.hpp" +#include "gemmi/it92.hpp" // IT92 x-ray form factors +#include "gemmi/model.hpp" + +// gemmi's gridding of a model - the density of DensityCalculator::put_model_density_on_grid() and the +// bulk-solvent mask of SolventMasker(AtomicRadiiSet::Refmac).put_mask_on_grid() - spread over +// `nthreads` threads of ParallelFor's pool. Every grid point receives the same operations in the same +// order as in gemmi's serial loops, so the grids are gemmi's, bit for bit, on any number of threads. + +// The points gemmi's Grid::symmetrize() reduces the orbits of this grid's size and space group into +// (empty in P1). Both functions below take them, for a grid of the same size. +std::vector OrbitLeaders(const gemmi::Grid &grid, size_t nthreads); + +// `dc` set up as for put_model_density_on_grid(): d_min, rate, blur, the grid's cell and group. +void PutModelDensityOnGrid(gemmi::DensityCalculator, float> &dc, const gemmi::Model &model, + const std::vector &orbit_leaders, size_t nthreads); + +// `grid` sized, with its cell and group set. +void PutMaskOnGrid(gemmi::Grid &grid, const gemmi::Model &model, const std::vector &orbit_leaders, + size_t nthreads); diff --git a/rugnux/RigidBodyRefine.cpp b/rugnux/RigidBodyRefine.cpp index e60339c57..148230039 100644 --- a/rugnux/RigidBodyRefine.cpp +++ b/rugnux/RigidBodyRefine.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -16,13 +17,12 @@ #include "gemmi/dencalc.hpp" // DensityCalculator #include "gemmi/it92.hpp" // IT92 x-ray form factors #include "gemmi/scaling.hpp" // Scaling (bulk solvent + anisotropic B) -#include "gemmi/solmask.hpp" // SolventMasker #include "ModelFFT.h" // MapToFPhi +#include "ModelGrid.h" // PutModelDensityOnGrid, PutMaskOnGrid #include "ModelScaling.h" // FitModelScale #include "../common/JFJochMath.h" // PI #include "../common/Logger.h" -#include "../common/ParallelFor.h" namespace { @@ -69,6 +69,18 @@ struct Placement { } }; +// How finely a zone's maps are sampled: DensityCalculator's rate, for a spacing of d_min / (2 * rate). +constexpr double GRID_RATE = 1.5; + +// An empty grid of the zone's size, the size DensityCalculator gives its own at GRID_RATE. +gemmi::Grid ZoneGrid(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min) { + gemmi::Grid grid; + grid.unit_cell = cell; + grid.spacegroup = &sg; + grid.set_size_from_spacing(d_min / (2 * GRID_RATE), gemmi::GridSizeRounding::Up); + return grid; +} + // One target evaluation: place the model, recompute Fcalc and the bulk-solvent mask to the zone's // resolution, re-fit the scale, and hand back the amplitude residuals. class Evaluator { @@ -77,10 +89,13 @@ public: const std::vector &base, const Placement &placement, size_t nthreads) : model_(model), cell_(cell), sg_(sg), base_(base), placement_(placement), nthreads_(nthreads) {} - // The zone's observations, and the scale the residuals are expressed in. - void SetZone(const gemmi::AsuData> &fobs, double d_min) { + // The zone's observations, and the scale the residuals are expressed in. `orbit_leaders` are those + // of the zone's grid, and must outlive the zone. + void SetZone(const gemmi::AsuData> &fobs, double d_min, + const std::vector &orbit_leaders) { fobs_ = fobs; d_min_ = d_min; + orbit_leaders_ = &orbit_leaders; double sum = 0; for (const auto &hv : fobs_.v) sum += hv.value.value; @@ -108,20 +123,16 @@ public: gemmi::DensityCalculator dc; dc.d_min = d_min_; - dc.rate = 1.5; + dc.rate = GRID_RATE; dc.grid.unit_cell = cell_; dc.grid.spacegroup = &sg_; dc.set_refmac_compatible_blur(model_); - dc.put_model_density_on_grid(model_); + PutModelDensityOnGrid(dc, model_, *orbit_leaders_, nthreads_); gemmi::AsuData> fcalc = MapToFPhi(dc.grid).prepare_asu_data(dc.d_min, dc.blur, false, false, false); - gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac); - gemmi::Grid mask_grid; - mask_grid.unit_cell = cell_; - mask_grid.spacegroup = &sg_; - mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up); - masker.put_mask_on_grid(mask_grid, model_); + gemmi::Grid mask_grid = ZoneGrid(cell_, sg_, d_min_); + PutMaskOnGrid(mask_grid, model_, *orbit_leaders_, nthreads_); gemmi::AsuData> fmask = MapToFPhi(mask_grid).prepare_asu_data(dc.d_min, 0); if (fmask.size() != fcalc.size()) @@ -187,6 +198,7 @@ private: Placement placement_; gemmi::AsuData> fobs_; double d_min_ = 0; + const std::vector *orbit_leaders_ = nullptr; double f_mean_ = 1; bool solvent_fitted_ = false; size_t nthreads_ = 1; @@ -204,6 +216,12 @@ private: // pair the central evaluation fitted - which is the pair the serial loop's shifted evaluations used, // since the zone's solvent is fitted by the first evaluation of the zone and every Evaluate starts // with the central one. Each column's arithmetic is the serial one, so the Jacobian is too. +// +// The columns run on threads of their own rather than on ParallelFor's pool, for the reason the null's +// replicates do (ModelValidation.cpp): an evaluation spreads its gridding over the pool, which it +// cannot do from a pool worker, where a parallel pass runs inline. On the pool, the six columns would +// each grid on one core and leave the rest of the machine idle for most of the fit. A run given one +// thread evaluates them one after the other, on its own. class RigidBodyCost : public ceres::CostFunction { public: RigidBodyCost(Evaluator &ev, std::vector &columns, size_t nthreads) @@ -229,13 +247,18 @@ public: const double step = ev_.JacobianStep(); std::vector> shifted(6, std::vector(n)); std::array ok{}; - ParallelFor(6, nthreads_, [&](int j) { - double q[6]; - std::copy(parameters[0], parameters[0] + 6, q); - q[j] += step; - columns_[j].UseSolvent(ev_.k_sol, ev_.b_sol); - ok[j] = columns_[j].Residuals(q, shifted[j].data()); - }); + const std::launch policy = nthreads_ > 1 ? std::launch::async : std::launch::deferred; + std::vector> running; + for (int j = 0; j < 6; j++) + running.push_back(std::async(policy, [&, j] { + double q[6]; + std::copy(parameters[0], parameters[0] + 6, q); + q[j] += step; + columns_[j].UseSolvent(ev_.k_sol, ev_.b_sol); + ok[j] = columns_[j].Residuals(q, shifted[j].data()); + })); + for (std::future &f : running) + f.get(); // Counted as the serial loop counted them: it stopped at the first column that failed. for (int j = 0; j < 6; j++) { ++ev_.evaluations; @@ -349,7 +372,7 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, std::vector columns; columns.reserve(6); for (int j = 0; j < 6; j++) - columns.emplace_back(column_models[j], cell, sg, base, placement, 1); + columns.emplace_back(column_models[j], cell, sg, base, placement, nthreads); double q[6] = {0, 0, 0, 0, 0, 0}; bool any_zone_solved = false; for (double zone : ladder) { @@ -362,9 +385,10 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model, if (zone_obs.v.size() < 50) continue; result.zones.push_back(zone); // the ladder WALKED, which a thin zone drops out of - ev.SetZone(zone_obs, zone); + const std::vector orbit_leaders = OrbitLeaders(ZoneGrid(cell, sg, zone), nthreads); + ev.SetZone(zone_obs, zone, orbit_leaders); for (Evaluator &c : columns) - c.SetZone(zone_obs, zone); + c.SetZone(zone_obs, zone, orbit_leaders); ceres::Problem problem; problem.AddResidualBlock(new RigidBodyCost(ev, columns, nthreads), nullptr, q); diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp index 69ead8d31..60703b4b9 100644 --- a/tests/ModelValidationTest.cpp +++ b/tests/ModelValidationTest.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include #include @@ -13,9 +14,11 @@ #include #include +#include #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 &a, const gemmi::Grid &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, 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 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 leaders = OrbitLeaders(gemmi_mask, nthreads); + + gemmi::DensityCalculator, 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 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)); + } + } + } +}