Merge branch 'model-grid-parallel' into rc173 (rigid body: density gridding, symmetrize and mask parallel and bit-identical)

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01D1G8gJVAy6gp1K5Dz3NE5C

docs/CHANGELOG.md

@@ -9,6 +9,7 @@
 * jfjoch_viewer's dark theme no longer leaves navy buttons, red warnings and chart guide lines at their light-theme colours.
 * `rugnux --model` draws its null's random placements away from every orientation equivalent to the model's under the space group or a twin law, so a random placement can no longer be refined onto the model's own solution and make the crystal's own model read as not fitting.
 * 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
 

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

rugnux/ModelGrid.cpp

@@ -0,0 +1,245 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include "ModelGrid.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <utility>
+
+#include "gemmi/solmask.hpp"   // SolventMasker
+
+#include "../common/ParallelFor.h"
+
+namespace {
+
+using Table = gemmi::IT92<float>;
+
+// The box of grid points gemmi's Grid::use_points_in_box<true>() 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<true>(fpos, radius, ..., false), sized and clamped as gemmi does it.
+AtomBox MakeAtomBox(const gemmi::Grid<float> &grid, const gemmi::Fractional &fpos, double radius) {
+    AtomBox box;
+    box.fpos = fpos;
+    box.radius = radius;
+    box.du = std::min(static_cast<int>(std::ceil(radius / grid.spacing[0])), grid.nu - 1);
+    box.dv = std::min(static_cast<int>(std::ceil(radius / grid.spacing[1])), grid.nv - 1);
+    box.dw = std::min(static_cast<int>(std::ceil(radius / grid.spacing[2])), grid.nw - 1);
+    return box;
+}
+
+// gemmi's use_points_in_box<true>() 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 <typename FuncFor>
+void UseAtomBoxesByPlane(gemmi::Grid<float> &grid, const std::vector<AtomBox> &boxes, size_t nthreads,
+                         FuncFor func_for) {
+    const int nu = grid.nu, nv = grid.nv, nw = grid.nw;
+    std::vector<std::vector<int>> plane_atoms(nw);
+    for (int i = 0; i < static_cast<int>(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<int> &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 <typename Func>
+void SymmetrizeOrbits(gemmi::Grid<float> &grid, const std::vector<size_t> &leaders, size_t nthreads,
+                      Func func) {
+    const std::vector<gemmi::GridOp> ops = grid.get_scaled_ops_except_id();
+    const size_t nu = grid.nu, nv = grid.nv;
+    ParallelChunks(static_cast<int>(leaders.size()), nthreads, [&](int lo, int hi) {
+        std::vector<size_t> mates(ops.size());
+        for (int i = lo; i < hi; i++) {
+            const size_t idx = leaders[i];
+            const int u = static_cast<int>(idx % nu), v = static_cast<int>(idx / nu % nv),
+                      w = static_cast<int>(idx / (nu * nv));
+            for (size_t k = 0; k < ops.size(); ++k) {
+                const std::array<int, 3> 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<const gemmi::Atom *> ModelAtoms(const gemmi::Model &model) {
+    std::vector<const gemmi::Atom *> 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<const Table::Coef &>().precalculate_density_iso(CReal(), CReal()));
+    using AnisoSum = decltype(std::declval<const Table::Coef &>().precalculate_density_aniso_b(
+        gemmi::SMat33<CReal>(), 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<size_t> OrbitLeaders(const gemmi::Grid<float> &grid, size_t nthreads) {
+    const std::vector<gemmi::GridOp> ops = grid.get_scaled_ops_except_id();
+    if (ops.empty())
+        return {};   // P1: nothing to symmetrize
+    std::vector<std::vector<size_t>> 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<int, 3> 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<size_t> leaders;
+    for (const std::vector<size_t> &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<Table, float> &dc, const gemmi::Model &model,
+                           const std::vector<size_t> &orbit_leaders, size_t nthreads) {
+    using CReal = AtomDensity::CReal;
+    dc.initialize_grid();
+    const std::vector<const gemmi::Atom *> atoms = ModelAtoms(model);
+    const int n = static_cast<int>(atoms.size());
+    std::vector<AtomBox> boxes(n);
+    std::vector<AtomDensity> 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<CReal>(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<float> &grid, const gemmi::Model &model, const std::vector<size_t> &orbit_leaders,
+                   size_t nthreads) {
+    const gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
+    masker.clear(grid);
+    std::vector<AtomBox> 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);
+}

rugnux/ModelGrid.h

@@ -0,0 +1,29 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#pragma once
+
+#include <cstddef>
+#include <vector>
+
+#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<size_t> OrbitLeaders(const gemmi::Grid<float> &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<gemmi::IT92<float>, float> &dc, const gemmi::Model &model,
+                           const std::vector<size_t> &orbit_leaders, size_t nthreads);
+
+// `grid` sized, with its cell and group set.
+void PutMaskOnGrid(gemmi::Grid<float> &grid, const gemmi::Model &model, const std::vector<size_t> &orbit_leaders,
+                   size_t nthreads);

rugnux/RigidBodyRefine.cpp

@@ -8,6 +8,7 @@
 #include <chrono>
 #include <cmath>
 #include <complex>
+#include <future>
 #include <vector>
 
 #include <ceres/ceres.h>
@@ -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<float> ZoneGrid(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min) {
+    gemmi::Grid<float> 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<gemmi::Position> &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<gemmi::ValueSigma<float>> &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<gemmi::ValueSigma<float>> &fobs, double d_min,
+                 const std::vector<size_t> &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<Table, float> 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<std::complex<float>> fcalc =
             MapToFPhi(dc.grid).prepare_asu_data(dc.d_min, dc.blur, false, false, false);
 
-        gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
-        gemmi::Grid<float> 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<float> mask_grid = ZoneGrid(cell_, sg_, d_min_);
+        PutMaskOnGrid(mask_grid, model_, *orbit_leaders_, nthreads_);
         gemmi::AsuData<std::complex<float>> 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<gemmi::ValueSigma<float>> fobs_;
     double d_min_ = 0;
+    const std::vector<size_t> *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<Evaluator> &columns, size_t nthreads)
@@ -229,13 +247,18 @@ public:
             const double step = ev_.JacobianStep();
             std::vector<std::vector<double>> shifted(6, std::vector<double>(n));
             std::array<bool, 6> 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<std::future<void>> 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<void> &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;
@@ -367,7 +390,7 @@ RigidBodyRefineResult RefineRigidBody(gemmi::Model &model,
     std::vector<Evaluator> 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) {
@@ -380,9 +403,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<size_t> 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);

tests/ModelValidationTest.cpp

@@ -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));
+            }
+        }
+    }
+}