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