Merge ribo-gap into refine-merge: the bulk-solvent mask grid never coarser than 0.75 A, in refinement and validation

Resolved: the ModelGrid.h include in ModelStructureFactorsGPU.cpp (AnisoRadiusB, SolventMaskGrid).
ModelReferenceIntensities keeps its mask on the density's grid: it is a reference at fixed solvent
constants, not a scored mask, and RigidBodyGPU_MatchesCPU and ModelScaleGPU_MatchesGemmiOnAModelsPoints
make their amplitudes from it - with the finer mask those amplitudes moved enough to land the anisotropic
test cluster on a near-tie of gemmi's scale fit (residuals 2.8e-4 against 2e-4, b_star 9.5e-7 against
1.6e-7). ModelStructureFactorsGPU_MatchesCPU at 1.5 A in F 41 3 2 (F_calc 1.28e-4 of the mean |F|
against 1e-4) fails on 89b37f9bc already and is left as it is.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

docs/CPU_DATA_ANALYSIS_DECISIONS.md

@@ -132,7 +132,7 @@ The model electron density is sampled on a grid (IT92 X-ray form factors, with a
 
 ### 14.2 Bulk solvent and scaling
 
-A flat bulk-solvent mask around the model is transformed to $F_\mathrm{mask}$, and the model is scaled to the observed amplitudes by an overall least-squares fit of a scale $k$, an anisotropic $B$, and the flat-solvent parameters $k_\mathrm{sol}, B_\mathrm{sol}$:
+A flat bulk-solvent mask around the model (the REFMAC radii, on the density's grid, or on one of its own at 0.75 Å where that grid is too coarse for the mask's 0.8 Å shrink to act) is transformed to $F_\mathrm{mask}$, and the model is scaled to the observed amplitudes by an overall least-squares fit of a scale $k$, an anisotropic $B$, and the flat-solvent parameters $k_\mathrm{sol}, B_\mathrm{sol}$:
 
 $
 F_\mathrm{model} = k\,e^{-\mathbf{h}^\top \mathbf{B}\,\mathbf{h}/4}\left(F_\mathrm{calc} + k_\mathrm{sol}\,e^{-B_\mathrm{sol}\,s^2}\,F_\mathrm{mask}\right),\quad s^2 = 1/4d^2.

docs/CPU_DATA_ANALYSIS_REFINEMENT.md

@@ -93,8 +93,8 @@ gradient gather, whatever the number of atoms. The derivative of a riding hydrog
 
 **Scale and solvent.** The flat bulk-solvent mask is drawn with the REFMAC radii (probe 1.0 Å, shrink 0.8 Å)
 on a grid of its own, with a spacing of $d_\mathrm{min}/3$ or 0.4 Å, whichever is finer (coarser only where
-a cell would need more than $2^{26}$ points, and then never coarser than the density grid, which a cell that
-large uses for its mask). On the atoms' density grid the shrink is lost at low
+a cell would need more than $2^{26}$ points, and then never coarser than 0.75 Å nor than the density grid,
+which a cell that large uses for its mask). On the atoms' density grid the shrink is lost at low
 resolution: at 3 Å its 1 Å step leaves no grid point within the shrink radius. $k_\mathrm{sol}$ and
 $B_\mathrm{sol}$ are searched over their physical range on the working reflections, and the overall scale
 and anisotropic B are fitted at each point of that search (§14.2) - at the points of its fine pass by

docs/GPU_DATA_ANALYSIS.md

@@ -166,7 +166,8 @@ the first card has no room for go on the others, each against its own budget.
   as a gather over bricks of the grid, bit-identical from run to run.
   [§14.1](CPU_DATA_ANALYSIS_DECISIONS.md#141-model-structure-factors)
 - **Bulk-solvent mask** (`ModelMaskGPU.cu`) - the solvent mask from every symmetry image of every
-  atom, with island removal and the shrink, as gemmi computes it.
+  atom, with island removal and the shrink, as gemmi computes it - on the density's grid, or on one of its
+  own at 0.75 A where the density's is too coarse for the 0.8 A shrink to act.
   [§14.2](CPU_DATA_ANALYSIS_DECISIONS.md#142-bulk-solvent-and-scaling)
 - **Structure factors and maps** (`ModelStructureFactorsGPU.cu`) - F_calc and F_mask from the two grids
   by cuFFT, composed over the group's operators and read off at the data's reflections; and the
@@ -269,7 +270,8 @@ and both gradients come back into pinned buffers. The minimiser (L-BFGS) stays o
 Once per macro-cycle ([§15.3](CPU_DATA_ANALYSIS_REFINEMENT.md#153-scale-bulk-solvent-and-the-likelihood-parameters)):
 
 - **F_calc and F_mask** - the forward path above without the target, and the bulk-solvent mask drawn on the card
-  (`ModelMaskGPU.cu`) on **its own fine grid** (d_min/3 or 0.4 A, whichever is finer, at most 2^26 points but never coarser than the density's grid),
+  (`ModelMaskGPU.cu`) on **its own fine grid** (d_min/3 or 0.4 A, whichever is finer, at most 2^26 points but never coarser than 0.75 A nor than the
+  density's grid),
   with a grid, spectrum and cuFFT plan of its own where that grid is not the density's: one block per
   (atom, symmetry image) sets the points within the atom's REFMAC radius, the islands are removed by a
   union-find labelling (integer `atomicCAS`, rooted at each component's lowest index, so the partition does not

image_analysis/structure_refinement/ModelGrid.cpp

@@ -501,6 +501,23 @@ void PutMaskOnGrid(gemmi::Grid<float> &grid, const gemmi::Model &model, const st
         masker.shrink(grid);
 }
 
+gemmi::Grid<float> SolventMaskGrid(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min,
+                                   double spacing) {
+    spacing = std::min({spacing, MASK_SPACING_SHRINK, d_min / 3.0});
+    gemmi::Grid<float> grid;
+    grid.unit_cell = cell;
+    grid.spacegroup = &sg;
+    // set_size_from_spacing()'s size, without allocating the values
+    const std::array<int, 3> n = gemmi::good_grid_size({{cell.a / spacing, cell.b / spacing, cell.c / spacing}},
+                                                       gemmi::GridSizeRounding::Up, &sg);
+    grid.nu = n[0];
+    grid.nv = n[1];
+    grid.nw = n[2];
+    grid.calculate_spacing();
+    grid.axis_order = gemmi::AxisOrder::XYZ;
+    return grid;
+}
+
 size_t ModelGridPoints(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg, double d_min, double rate) {
     const double spacing = d_min / (2 * rate);
     const std::array<int, 3> m = gemmi::good_grid_size({{cell.a / spacing, cell.b / spacing, cell.c / spacing}},

image_analysis/structure_refinement/ModelGrid.h

@@ -61,6 +61,17 @@ void PutModelDensityOnGrid(gemmi::DensityCalculator<gemmi::IT92<float>, float> &
 void PutMaskOnGrid(gemmi::Grid<float> &grid, const gemmi::Model &model, const std::vector<size_t> &orbit_leaders,
                    size_t nthreads);
 
+// The step a bulk-solvent mask's grid is held below. The mask's shrink (0.8 A, the Refmac radii) reaches no grid point
+// on a coarser grid and is skipped, and the macromolecule keeps the whole 1.0 A probe: some 5-10% of the cell's solvent
+// lost, 0.015-0.03 of R-free on deposited models. Below the shrink radius the R-factors no longer depend on the step.
+constexpr double MASK_SPACING_SHRINK = 0.75;   // A
+
+// The grid of a bulk-solvent mask for structure factors to d_min (sizes, cell and group set, no values): at `spacing`,
+// but never coarser than MASK_SPACING_SHRINK nor than the density's own grid at rate 1.5 (d_min / 3), which it then
+// is - set_size_from_spacing()'s size for d_min / 3.
+gemmi::Grid<float> SolventMaskGrid(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min,
+                                   double spacing);
+
 // The number of points of the grid a model's structure factors to d_min are computed on at DensityCalculator's
 // `rate` (a spacing of d_min / (2 rate)) in this cell and space group: the grid's size, without making it.
 size_t ModelGridPoints(const gemmi::UnitCell &cell, const gemmi::SpaceGroup *sg, double d_min, double rate = 1.5);

image_analysis/structure_refinement/ModelStructureFactorsGPU.cpp

@@ -13,7 +13,7 @@
 #include "gemmi/fourier.hpp"   // get_size_for_hkl, get_f_phi_on_grid
 #include "gemmi/solmask.hpp"   // SolventMasker, refmac_radius_for_bulk_solvent
 
-#include "ModelGrid.h"         // AnisoRadiusB
+#include "ModelGrid.h"         // AnisoRadiusB, SolventMaskGrid
 #include "RigidBodyRefine.h"   // RigidBodyZoneGrid
 
 namespace {
@@ -123,6 +123,11 @@ ModelStructureFactorsGPU::ModelStructureFactorsGPU(int device, const gemmi::Mode
             setup_.grid.frac[3 * i + j] = cell.frac.mat[i][j];
         }
     setup_.volume = cell.volume;
+    const gemmi::Grid<float> mask_grid = SolventMaskGrid(cell, sg, d_min, d_min / 3.0);
+    setup_.mask_grid = setup_.grid;
+    setup_.mask_grid.nu = mask_grid.nu;
+    setup_.mask_grid.nv = mask_grid.nv;
+    setup_.mask_grid.nw = mask_grid.nw;
 
     const gemmi::GroupOps gops = sg.operations();
     setup_.den = gemmi::Op::DEN;
@@ -177,7 +182,7 @@ ModelStructureFactorsGPU::ModelStructureFactorsGPU(int device, const gemmi::Mode
     const int map_nw = 2 * (map_size_[2] / 2);
     setup_.map_points = static_cast<size_t>(map_size_[0]) * map_size_[1] * map_nw;
     setup_.map_complex_points = static_cast<size_t>(map_size_[0]) * map_size_[1] * (map_nw / 2 + 1);
-    fft_work_bytes_ = ModelStructureFactorsGPUEngine::FFTWorkBytes(device_, setup_.grid, map_size_[0], map_size_[1], map_nw);
+    fft_work_bytes_ = ModelStructureFactorsGPUEngine::FFTWorkBytes(device_, setup_, map_size_[0], map_size_[1], map_nw);
     device_bytes_ = ModelStructureFactorsGPUEngine::DeviceBytes(setup_, fft_work_bytes_);
 }
 

image_analysis/structure_refinement/ModelStructureFactorsGPU.cu

@@ -151,6 +151,10 @@ size_t ComplexPoints(const ModelDensityGrid &g) {
     return static_cast<size_t>(g.nu / 2 + 1) * g.nv * g.nw;
 }
 
+bool SameGrid(const ModelDensityGrid &a, const ModelDensityGrid &b) {
+    return a.nu == b.nu && a.nv == b.nv && a.nw == b.nw;
+}
+
 }  // namespace
 
 struct ModelStructureFactorsGPUEngineImpl {
@@ -161,7 +165,7 @@ struct ModelStructureFactorsGPUEngineImpl {
 
     CudaDevicePtr<float4> pos;
     CudaDevicePtr<ModelMaskAtom> mask_atoms;
-    CudaDevicePtr<float> grid;          // the density, then the mask; a map's real grid
+    CudaDevicePtr<float> grid;          // the density, then the mask (on a grid of its own size or not); a map's real grid
     CudaDevicePtr<float2> spectrum;     // their transforms; a map's coefficients, transformed in place
     CudaDevicePtr<char> fft_work;
     CudaDevicePtr<int> ops, row_hkl;
@@ -169,7 +173,7 @@ struct ModelStructureFactorsGPUEngineImpl {
     CudaDevicePtr<float2> f, fmask;
     std::unique_ptr<ModelDensityGPU> density;
     std::unique_ptr<ModelMaskGPU> mask;
-    CudaFFTPlan forward;
+    CudaFFTPlan forward, mask_forward;   // the mask's only where its grid is not the density's
 
     // On `dev`, which the caller has made current.
     ModelStructureFactorsGPUEngineImpl(int dev, const ModelStructureFactorsGPUSetup &s, size_t work)
@@ -178,8 +182,9 @@ struct ModelStructureFactorsGPUEngineImpl {
         const size_t na = std::max<size_t>(setup.max_atoms, 1), nr = std::max<size_t>(setup.rows.size(), 1);
         pos = CudaDevicePtr<float4>(na, sync);
         mask_atoms = CudaDevicePtr<ModelMaskAtom>(na, sync);
-        grid = CudaDevicePtr<float>(std::max(GridPoints(setup.grid), setup.map_points), sync);
-        spectrum = CudaDevicePtr<float2>(std::max(ComplexPoints(setup.grid), setup.map_complex_points), sync);
+        grid = CudaDevicePtr<float>(std::max({GridPoints(setup.grid), GridPoints(setup.mask_grid), setup.map_points}), sync);
+        spectrum = CudaDevicePtr<float2>(
+            std::max({ComplexPoints(setup.grid), ComplexPoints(setup.mask_grid), setup.map_complex_points}), sync);
         fft_work = CudaDevicePtr<char>(std::max<size_t>(work_bytes, 1), sync);
         ops = CudaDevicePtr<int>(std::max<size_t>(12 * setup.sym_ops.size(), 1), sync);
         row_hkl = CudaDevicePtr<int>(3 * nr, sync);
@@ -188,13 +193,13 @@ struct ModelStructureFactorsGPUEngineImpl {
         fmask = CudaDevicePtr<float2>(nr, sync);
         density = std::make_unique<ModelDensityGPU>(stream, setup.max_atoms, setup.max_pairs,
                                                     ModelDensityGPU::Bricks(setup.grid.nu, setup.grid.nv, setup.grid.nw));
-        mask = std::make_unique<ModelMaskGPU>(stream, GridPoints(setup.grid));
+        mask = std::make_unique<ModelMaskGPU>(stream, GridPoints(setup.mask_grid));
 
         ModelMaskGrid mask_grid;
-        mask_grid.nu = setup.grid.nu;
-        mask_grid.nv = setup.grid.nv;
-        mask_grid.nw = setup.grid.nw;
-        std::copy(setup.grid.orth, setup.grid.orth + 9, mask_grid.orth);
+        mask_grid.nu = setup.mask_grid.nu;
+        mask_grid.nv = setup.mask_grid.nv;
+        mask_grid.nw = setup.mask_grid.nw;
+        std::copy(setup.mask_grid.orth, setup.mask_grid.orth + 9, mask_grid.orth);
         mask_grid.volume = setup.volume;
         mask->SetGrid(mask_grid, setup.mask_ops);
 
@@ -207,7 +212,10 @@ struct ModelStructureFactorsGPUEngineImpl {
         up(row_scale, setup.row_scale.data(), setup.row_scale.size() * sizeof(double));
 
         forward = CudaFFTPlan(ForwardShape(setup.grid), stream, fft_work);
-        if (forward.WorkBytes() > work_bytes)
+        if (!SameGrid(setup.mask_grid, setup.grid))
+            mask_forward = CudaFFTPlan(ForwardShape(setup.mask_grid), stream, fft_work);
+        if (forward.WorkBytes() > work_bytes ||
+            (!SameGrid(setup.mask_grid, setup.grid) && mask_forward.WorkBytes() > work_bytes))
             throw JFJochException(JFJochExceptionCategory::GPUCUDAError, "structure factors: cuFFT work area over the reserve");
         cuda_err(cudaStreamSynchronize(stream));
     }
@@ -216,8 +224,8 @@ struct ModelStructureFactorsGPUEngineImpl {
         cudaStreamSynchronize(stream);
     }
 
-    float Norm() const {
-        return static_cast<float>(setup.volume / static_cast<double>(GridPoints(setup.grid)));
+    float Norm(const ModelDensityGrid &g) const {
+        return static_cast<float>(setup.volume / static_cast<double>(GridPoints(g)));
     }
 };
 
@@ -228,22 +236,25 @@ size_t ModelStructureFactorsGPUEngine::TotalMemory(int device) {
     return total;
 }
 
-size_t ModelStructureFactorsGPUEngine::FFTWorkBytes(int device, const ModelDensityGrid &grid, int map_nu, int map_nv,
-                                                    int map_nw) {
+size_t ModelStructureFactorsGPUEngine::FFTWorkBytes(int device, const ModelStructureFactorsGPUSetup &setup, int map_nu,
+                                                    int map_nv, int map_nw) {
     DeviceGuard guard(device);
-    return std::max(CudaFFTPlan::WorkBytes(ForwardShape(grid)), CudaFFTPlan::WorkBytes(MapShape(map_nu, map_nv, map_nw)));
+    size_t bytes = std::max(CudaFFTPlan::WorkBytes(ForwardShape(setup.grid)),
+                            CudaFFTPlan::WorkBytes(MapShape(map_nu, map_nv, map_nw)));
+    if (!SameGrid(setup.mask_grid, setup.grid))
+        bytes = std::max(bytes, CudaFFTPlan::WorkBytes(ForwardShape(setup.mask_grid)));
+    return bytes;
 }
 
 size_t ModelStructureFactorsGPUEngine::DeviceBytes(const ModelStructureFactorsGPUSetup &s, size_t fft_work_bytes) {
     const size_t na = s.max_atoms + 1, nr = s.rows.size() + 1;
-    const size_t points = GridPoints(s.grid);
     size_t bytes = na * (sizeof(float4) + sizeof(ModelMaskAtom));
-    bytes += std::max(points, s.map_points) * sizeof(float);
-    bytes += std::max(ComplexPoints(s.grid), s.map_complex_points) * sizeof(float2);
+    bytes += std::max({GridPoints(s.grid), GridPoints(s.mask_grid), s.map_points}) * sizeof(float);
+    bytes += std::max({ComplexPoints(s.grid), ComplexPoints(s.mask_grid), s.map_complex_points}) * sizeof(float2);
     bytes += fft_work_bytes;
     bytes += 12 * s.sym_ops.size() * sizeof(int) + nr * (3 * sizeof(int) + sizeof(double) + 2 * sizeof(float2));
     bytes += ModelDensityGPU::DeviceBytes(s.max_atoms, s.max_pairs, ModelDensityGPU::Bricks(s.grid.nu, s.grid.nv, s.grid.nw));
-    bytes += ModelMaskGPU::DeviceBytes(points);
+    bytes += ModelMaskGPU::DeviceBytes(GridPoints(s.mask_grid));
     return bytes;
 }
 
@@ -282,17 +293,19 @@ void ModelStructureFactorsGPUEngine::Compute(const std::vector<ModelDensityAtom>
     e.density->Compute(e.pos, e.grid);
     cufft_err(e.forward.R2C(e.grid, reinterpret_cast<cufftComplex *>(e.spectrum.get())));
     if (rows > 0) {
-        compose_kernel<<<blocks(rows), THREADS, 0, e.stream>>>(e.spectrum, e.Norm(), g.nu, g.nv, g.nw, e.ops,
+        compose_kernel<<<blocks(rows), THREADS, 0, e.stream>>>(e.spectrum, e.Norm(g), g.nu, g.nv, g.nw, e.ops,
                                                                static_cast<int>(e.setup.sym_ops.size()), e.setup.den,
                                                                e.row_hkl, e.row_scale, rows, e.f);
         cuda_err(cudaGetLastError());
     }
 
+    const ModelDensityGrid &mg = e.setup.mask_grid;
     e.mask->Compute(e.mask_atoms, static_cast<int>(mask_atoms.size()), e.grid);
-    cufft_err(e.forward.R2C(e.grid, reinterpret_cast<cufftComplex *>(e.spectrum.get())));
+    CudaFFTPlan &mask_forward = SameGrid(mg, g) ? e.forward : e.mask_forward;
+    cufft_err(mask_forward.R2C(e.grid, reinterpret_cast<cufftComplex *>(e.spectrum.get())));
     if (rows > 0) {
-        fmask_kernel<<<blocks(rows), THREADS, 0, e.stream>>>(e.spectrum, e.Norm(), e.row_hkl, rows, g.nu, g.nv, g.nw,
-                                                             e.fmask);
+        fmask_kernel<<<blocks(rows), THREADS, 0, e.stream>>>(e.spectrum, e.Norm(mg), e.row_hkl, rows, mg.nu, mg.nv,
+                                                             mg.nw, e.fmask);
         cuda_err(cudaGetLastError());
     }
 

image_analysis/structure_refinement/ModelStructureFactorsGPUEngine.h

@@ -19,6 +19,7 @@
 // Everything an engine needs that does not depend on the model's coordinates.
 struct ModelStructureFactorsGPUSetup {
     ModelDensityGrid grid;                       // the structure-factor grid, u fastest
+    ModelDensityGrid mask_grid;                  // the bulk-solvent mask's (SolventMaskGrid, ModelGrid.h): `grid`, or finer
     double volume = 0;
     // The group's operators as gemmi holds them: rot (row-major) and tran, in units of 1 / den.
     std::vector<std::array<int, 12>> sym_ops;
@@ -39,8 +40,9 @@ public:
 
     // The device's total memory.
     static size_t TotalMemory(int device);
-    // The cuFFT work area the structure factors' r2c, and a map of (nu, nv, nw)'s c2r, take on `device`.
-    static size_t FFTWorkBytes(int device, const ModelDensityGrid &grid, int map_nu, int map_nv, int map_nw);
+    // The cuFFT work area the structure factors' and the mask's r2c, and a map of (nu, nv, nw)'s c2r, take on `device`.
+    static size_t FFTWorkBytes(int device, const ModelStructureFactorsGPUSetup &setup, int map_nu, int map_nv,
+                               int map_nw);
     // The bytes an engine for `setup` reserves on the device, given its cuFFT work area.
     static size_t DeviceBytes(const ModelStructureFactorsGPUSetup &setup, size_t fft_work_bytes);
 

image_analysis/structure_refinement/ModelValidation.cpp

@@ -896,18 +896,24 @@ ModelValidationResult Validate(const std::vector<MergedReflection> &merged,
         // The mask is made on a thread of its own, beside the density: both only read the model, each
         // fills a grid of its own, and between them they are the whole cost of this step. Both are
         // gemmi's gridding spread over the threads (ModelGrid.h), the same grids bit for bit at any
-        // thread count; the two grids are the same size, so they share the symmetry's orbits.
+        // thread count. The mask is on the density's grid, and shares its symmetry's orbits, but where that
+        // grid is too coarse for the mask's shrink (SolventMaskGrid).
         // Refmac radii give a slightly lower R than the Cctbx set on our test cases, at no cost.
-        gemmi::Grid<float> mask_grid;
-        mask_grid.unit_cell = dc.grid.unit_cell;
-        mask_grid.spacegroup = dc.grid.spacegroup;
-        mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
-        const std::vector<size_t> orbit_leaders = OrbitLeaders(mask_grid, nthreads);
+        gemmi::Grid<float> mask_grid =
+            SolventMaskGrid(dc.grid.unit_cell, *dc.grid.spacegroup, to_d, dc.requested_grid_spacing());
+        mask_grid.data.resize(mask_grid.point_count());
+        gemmi::Grid<float> sized;
+        sized.unit_cell = dc.grid.unit_cell;
+        sized.spacegroup = dc.grid.spacegroup;
+        sized.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
+        const bool same = sized.nu == mask_grid.nu && sized.nv == mask_grid.nv && sized.nw == mask_grid.nw;
+        const std::vector<size_t> mask_leaders = OrbitLeaders(mask_grid, nthreads);
+        const std::vector<size_t> density_leaders = same ? mask_leaders : OrbitLeaders(sized, nthreads);
         std::future<gemmi::AsuData<std::complex<float>>> fmask = std::async(std::launch::async, [&] {
-            PutMaskOnGrid(mask_grid, ms.st.models[0], orbit_leaders, nthreads);
+            PutMaskOnGrid(mask_grid, ms.st.models[0], mask_leaders, nthreads);
             return MapToFPhi(mask_grid).prepare_asu_data(to_d, 0);
         });
-        PutModelDensityOnGrid(dc, ms.st.models[0], orbit_leaders, nthreads);
+        PutModelDensityOnGrid(dc, ms.st.models[0], density_leaders, nthreads);
         ms.fcalc = MapToFPhi(dc.grid)
                        .prepare_asu_data(dc.d_min, dc.blur, false, false, false);
         ms.fmask = fmask.get();
@@ -2174,14 +2180,18 @@ ModelValidationResult ReadModelAnomalousMap(const std::vector<MergedReflection>
         dc.rate = 1.5;
         dc.set_grid_cell_and_spacegroup(st);
         dc.set_refmac_compatible_blur(st.models[0]);
-        gemmi::Grid<float> mask_grid;
-        mask_grid.unit_cell = dc.grid.unit_cell;
-        mask_grid.spacegroup = dc.grid.spacegroup;
-        mask_grid.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
-        const std::vector<size_t> orbit_leaders = OrbitLeaders(mask_grid, nthreads);
-        PutMaskOnGrid(mask_grid, st.models[0], orbit_leaders, nthreads);
+        gemmi::Grid<float> mask_grid =
+            SolventMaskGrid(dc.grid.unit_cell, *dc.grid.spacegroup, d_min, dc.requested_grid_spacing());
+        mask_grid.data.resize(mask_grid.point_count());
+        gemmi::Grid<float> sized;
+        sized.unit_cell = dc.grid.unit_cell;
+        sized.spacegroup = dc.grid.spacegroup;
+        sized.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
+        const bool same = sized.nu == mask_grid.nu && sized.nv == mask_grid.nv && sized.nw == mask_grid.nw;
+        const std::vector<size_t> mask_leaders = OrbitLeaders(mask_grid, nthreads);
+        PutMaskOnGrid(mask_grid, st.models[0], mask_leaders, nthreads);
         fmask = MapToFPhi(mask_grid).prepare_asu_data(d_min, 0);
-        PutModelDensityOnGrid(dc, st.models[0], orbit_leaders, nthreads);
+        PutModelDensityOnGrid(dc, st.models[0], same ? mask_leaders : OrbitLeaders(sized, nthreads), nthreads);
         fcalc = MapToFPhi(dc.grid).prepare_asu_data(dc.d_min, dc.blur, false, false, false);
     }
 
@@ -2287,6 +2297,8 @@ std::vector<MergedReflection> ModelReferenceIntensities(const RugnuxModel &model
     // and a correlation on raw intensities would be led by them.
     constexpr double K_SOL = 0.35;
     constexpr double B_SOL = 46.0;
+    // On the density's grid, unlike the scored masks (SolventMaskGrid): a reference at fixed solvent constants,
+    // which the tests' amplitudes are made from as well.
     gemmi::SolventMasker masker(gemmi::AtomicRadiiSet::Refmac);
     gemmi::Grid<float> mask_grid;
     mask_grid.unit_cell = dc.grid.unit_cell;

image_analysis/structure_refinement/RefinementXray.cpp

@@ -444,29 +444,11 @@ BulkSolventMask::BulkSolventMask(const gemmi::UnitCell &cell, const gemmi::Space
 
 double BulkSolventMask::MaskSpacing(const gemmi::UnitCell &cell, double d_min) {
     const double spacing = std::max(std::min(d_min / 3.0, MASK_SPACING_MAX), std::cbrt(cell.volume / MASK_POINTS_MAX));
-    return std::min(spacing, d_min / 3.0);
+    return std::min({spacing, MASK_SPACING_SHRINK, d_min / 3.0});
 }
 
 gemmi::Grid<float> BulkSolventMask::MaskGrid(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min) {
-    // the cap would coarsen it beyond the atoms' grid: that grid, whose points the density holds anyway
-    if (std::cbrt(cell.volume / MASK_POINTS_MAX) > d_min / 3.0) {
-        gemmi::Grid<float> grid = RigidBodyZoneGrid(cell, sg, d_min);
-        std::vector<float>().swap(grid.data);
-        return grid;
-    }
-    gemmi::Grid<float> grid;
-    grid.unit_cell = cell;
-    grid.spacegroup = &sg;
-    const double spacing = MaskSpacing(cell, d_min);
-    // set_size_from_spacing()'s size, without allocating the values
-    const std::array<int, 3> n = gemmi::good_grid_size({{cell.a / spacing, cell.b / spacing, cell.c / spacing}},
-                                                       gemmi::GridSizeRounding::Up, &sg);
-    grid.nu = n[0];
-    grid.nv = n[1];
-    grid.nw = n[2];
-    grid.calculate_spacing();
-    grid.axis_order = gemmi::AxisOrder::XYZ;
-    return grid;
+    return SolventMaskGrid(cell, sg, d_min, MaskSpacing(cell, d_min));
 }
 
 void BulkSolventMask::Fmask(const gemmi::Model &model, const std::vector<gemmi::Miller> &hkl,

image_analysis/structure_refinement/RefinementXray.h

@@ -90,11 +90,10 @@ struct XrayRowModel {
 // probe and shrink radii (1.0 and 0.8 A), which a grid at d_min / 3 cannot resolve at low resolution - at 3 A
 // its 1 A step leaves no grid point within the shrink radius and the shrink is skipped. The step is therefore
 // d_min / 3 or MASK_SPACING_MAX, whichever is finer - where the mask grid stays within MASK_POINTS_MAX points; a
-// cell too large for that gets the spacing that fills it, but never one coarser than the atoms' density grid
-// (RigidBodyZoneGrid, d_min / 3), which it then is: the density's points are on the card anyway, and the mask's
-// own grid only ever saved memory. The R-factors of low-resolution models settle at a step of 0.4 A (within 0.001
-// of 0.3 A) and move by up to 0.005 between 0.4 and 0.8 A; a cell coarsened by the cap past the shrink radius has
-// no shrink at all, its macromolecule swollen by the whole probe, which costs far more.
+// cell too large for that gets the spacing that fills it, but never one coarser than MASK_SPACING_SHRINK, where the
+// shrink still acts, nor than the atoms' density grid (d_min / 3), which it then is (SolventMaskGrid, ModelGrid.h).
+// The R-factors of low-resolution models settle at a step of 0.4 A (within 0.001 of 0.3 A) and move by up to 0.005
+// between 0.4 and 0.8 A; past the shrink radius there is no shrink at all, which costs far more.
 // The flat solvent model of Fokine & Urzhumtsev (2002) Acta Cryst. D58, 1387-1392
 // The probe-and-shrink mask of Jiang & Brunger (1994) J. Mol. Biol. 243, 100-115
 class BulkSolventMask {
@@ -104,7 +103,7 @@ public:
     BulkSolventMask(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min, size_t nthreads);
     // The mask's grid (sizes and spacing set, values not allocated) - also the GPU engine's (RefinementGPU.cpp).
     static gemmi::Grid<float> MaskGrid(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min);
-    // The step that grid is asked for.
+    // The step that grid is asked for (before good_grid_size() rounds it finer).
     static double MaskSpacing(const gemmi::UnitCell &cell, double d_min);
     // Fmask (unscaled) of `model` at `hkl`.
     void Fmask(const gemmi::Model &model, const std::vector<gemmi::Miller> &hkl,

tests/ModelValidationTest.cpp

@@ -1159,6 +1159,50 @@ TEST_CASE("ModelValidation_RigidBodySameOnAnyNumberOfThreads", "[ModelValidation
 // gemmi's but for the last bits of its single-precision walk (a point at a radius' edge can also fall on
 // the other side, adding or losing a term at the cutoff). Both must be the same on any number of threads,
 // for isotropic and anisotropic atoms and for groups with and without centring.
+// The bulk-solvent mask's grid: the density's (d_min / 3) where that is fine enough for the mask's 0.8 A shrink,
+// and finer where it is not - there the shrink is skipped on the density's grid and the macromolecule keeps the
+// whole probe, so less of the cell is solvent; on the mask's own grid the solvent is what a fine grid finds.
+TEST_CASE("SolventMaskGrid_ShrinkActsAtLowResolution", "[ModelValidation]") {
+    const auto path = WriteTemp("solvent_mask_grid_test.pdb",
+                                ClusterPdb("CRYST1   40.000   50.000   60.000  90.00  90.00  90.00 P 21 21 21    4\n").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);
+    auto sized = [&](double spacing) {
+        gemmi::Grid<float> g;
+        g.unit_cell = st.cell;
+        g.spacegroup = sg;
+        g.set_size_from_spacing(spacing, gemmi::GridSizeRounding::Up);
+        return g;
+    };
+    auto solvent = [&](gemmi::Grid<float> g) {
+        g.data.resize(g.point_count());
+        PutMaskOnGrid(g, st.models[0], OrbitLeaders(g, 2), 2);
+        double sum = 0;
+        for (float v : g.data)
+            sum += v;
+        return sum / static_cast<double>(g.data.size());
+    };
+
+    // 1.8 A: the density's grid, 0.6 A
+    const gemmi::Grid<float> fine = SolventMaskGrid(st.cell, *sg, 1.8, 1.8 / 3.0);
+    const gemmi::Grid<float> density_fine = sized(1.8 / 3.0);
+    CHECK(fine.nu == density_fine.nu);
+    CHECK(fine.nv == density_fine.nv);
+    CHECK(fine.nw == density_fine.nw);
+    CHECK(fine.data.empty());
+
+    // 4 A: the density's grid is 1.33 A, the mask's at most MASK_SPACING_SHRINK
+    const gemmi::Grid<float> coarse = SolventMaskGrid(st.cell, *sg, 4.0, 4.0 / 3.0);
+    for (int i = 0; i < 3; i++)
+        CHECK(coarse.spacing[i] <= MASK_SPACING_SHRINK);
+    const double on_density = solvent(sized(4.0 / 3.0)), on_mask = solvent(coarse), on_fine = solvent(sized(0.4));
+    INFO("solvent fraction on the density's grid " << on_density << ", the mask's " << on_mask << ", 0.4 A " << on_fine);
+    CHECK(on_mask > on_density + 0.01);
+    CHECK(std::fabs(on_mask - on_fine) < 0.005);
+}
+
 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",
@@ -1829,7 +1873,8 @@ TEST_CASE("RigidBodyGPU_CallsLeaveTheCallerOnItsCard", "[ModelValidation][gpu]")
 
 // The model's structure factors on the GPU against the CPU path model validation takes - gemmi's density and
 // mask on the symmetrized grid, transformed by FFTW and read off by prepare_asu_data() - in the five groups,
-// at a resolution where the mask's shrink step does nothing and at one where it does. The same reflections
+// at a resolution where the mask has a grid of its own (SolventMaskGrid) and at one where it is on the
+// density's. The same reflections
 // in the same order; F_calc to rounding; F_mask to the few grid points at an atom's mask radius that float
 // and double distances put on different sides. The same atoms twice give the same bits.
 TEST_CASE("ModelStructureFactorsGPU_MatchesCPU", "[ModelValidation][gpu]") {
@@ -1843,10 +1888,8 @@ TEST_CASE("ModelStructureFactorsGPU_MatchesCPU", "[ModelValidation][gpu]") {
             auto dc = ZoneDensity(st, model, d_min);
             dc.put_model_density_on_grid(model);
             const auto ref_fc = MapToFPhi(dc.grid).prepare_asu_data(d_min, dc.blur, false, false, false);
-            gemmi::Grid<float> mask;
-            mask.unit_cell = st.cell;
-            mask.spacegroup = sg;
-            mask.set_size_from_spacing(dc.requested_grid_spacing(), gemmi::GridSizeRounding::Up);
+            gemmi::Grid<float> mask = SolventMaskGrid(st.cell, *sg, d_min, dc.requested_grid_spacing());
+            mask.data.resize(mask.point_count());
             gemmi::SolventMasker(gemmi::AtomicRadiiSet::Refmac).put_mask_on_grid(mask, model);
             const auto ref_fm = MapToFPhi(mask).prepare_asu_data(d_min, 0);