diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md
index 9364cd64b..c58a481d4 100644
--- a/docs/ACKNOWLEDGEMENT.md
+++ b/docs/ACKNOWLEDGEMENT.md
@@ -557,6 +557,11 @@ N. W. Moriarty, M. Mustyakimov, T. C. Terwilliger, A. Urzhumtsev, P. H. Zwart an
 352-367 [doi:10.1107/S0907444912001308](https://doi.org/10.1107/S0907444912001308). R-free itself:
 A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal
 structures" (1992), Nature 355, 472-475 [doi:10.1038/355472a0](https://doi.org/10.1038/355472a0).
+The model the refinement scores is scaled to the data with the overall scale and anisotropic B fitted by
+robust least squares - iteratively reweighted, with the pseudo-Huber weight of each reflection's residual -
+so that a few reflections the model cannot explain do not set the scale of the rest: P. J. Huber, "Robust
+estimation of a location parameter" (1964), Ann. Math. Statist. 35, 73-101
+[doi:10.1214/aoms/1177703732](https://doi.org/10.1214/aoms/1177703732).
 
 **Unmodelled density** — the blobs `--refine` lists are the connected regions of the difference map above
 a cutoff, found by flood fill as Coot's "unmodelled blobs" finds them, in GEMMI's implementation
diff --git a/image_analysis/structure_refinement/ModelScaling.cpp b/image_analysis/structure_refinement/ModelScaling.cpp
index 6a099189b..022a9fa73 100644
--- a/image_analysis/structure_refinement/ModelScaling.cpp
+++ b/image_analysis/structure_refinement/ModelScaling.cpp
@@ -25,6 +25,8 @@
 
 namespace {
 
+constexpr int ROBUST_ROUNDS = 4;   // reweightings of the robust fit (ModelScaleBox::robust)
+
 // gemmi::Scaling<float> as a grid point fits it: the same starting point, parameters, model values,
 // derivatives and R. The solvent pair is fixed there, so |Fcalc + k_sol exp(-b_sol s^2) Fmask| of a
 // reflection is the same at every evaluation of the solver; it is taken once per grid point here rather than
@@ -35,8 +37,9 @@ struct FixedSolventFit {
         gemmi::Miller hkl;
         float fobs;
         float fcalc_abs;   // std::abs(Scaling::get_fcalc(p))
+        double weight = 1.0;   // 1, but for the robust fit's reweighting (RobustReweight)
         double get_y() const { return fobs; }
-        double get_weight() const { return 1.0; }
+        double get_weight() const { return weight; }
     };
     const std::vector<gemmi::Vec6> &constraint_matrix;
     double k_overall = 1.0;
@@ -61,7 +64,7 @@ struct FixedSolventFit {
             // Scaling::get_fcalc() at (k_sol, b_sol)
             const std::complex<float> fcalc =
                 s.use_solvent ? p.fcmol + (float) (k_sol * std::exp(-b_sol * p.stol2)) * p.fmask : p.fcmol;
-            points[i] = {p.hkl, p.fobs, std::abs(fcalc)};
+            points[i] = {p.hkl, p.fobs, std::abs(fcalc), 1.0};
             if (p.fobs < 1 || p.fobs < p.sigma)  // skip weak reflections
                 continue;
             double fcalc_abs = points[i].fcalc_abs;
@@ -138,6 +141,27 @@ struct FixedSolventFit {
     }
 };
 
+// The robust fit's weights at the current parameters: each reflection's residual r = |Fo| - |F_model| weighted by
+// 1 / sqrt(1 + (r / c)^2), with c the median |r| of the unweighted fit. Iterated, the weighted least squares then
+// minimises sum c^2 2 (sqrt(1 + (r / c)^2) - 1) - quadratic for the bulk of the reflections and linear for the few
+// a model fits badly (a beam-stop shadow, an ice ring, a missing part of the model), which least squares would
+// otherwise let set the scale of all the rest.
+// Following Huber (1964) Ann. Math. Statist. 35, 73-101 for the M-estimator, as iteratively reweighted least squares
+void RobustReweight(FixedSolventFit &fit, double &c) {
+    std::vector<double> r(fit.points.size());
+    for (size_t i = 0; i < fit.points.size(); i++)
+        r[i] = std::fabs(fit.points[i].fobs - fit.compute_value(fit.points[i]));
+    if (c <= 0) {
+        std::vector<double> sorted = r;
+        std::nth_element(sorted.begin(), sorted.begin() + sorted.size() / 2, sorted.end());
+        c = sorted[sorted.size() / 2];
+    }
+    if (c <= 0)
+        return;
+    for (size_t i = 0; i < fit.points.size(); i++)
+        fit.points[i].weight = 1.0 / std::sqrt(1.0 + gemmi::sq(r[i] / c));
+}
+
 // R-factor of the current parameters over the fitted reflections. This is what the grid is
 // selected on, and it is the quantity the scale exists to make small.
 double RFactor(const FixedSolventFit &fit) {
@@ -171,16 +195,17 @@ double ComputeLMMatrices(const FixedSolventFit &fit, const double *kaniso, doubl
         const double k = kaniso ? kaniso[i] : std::exp(-0.25 * fit.b_star.r_u_r(gemmi::Vec3(p.hkl)));
         double y = fit.compute_value_and_derivatives<NA>(p, k, dy_da);
         double dy_sig = p.get_y() - y;
+        const double w = p.get_weight();
         for (int j = 0; j != NA; ++j) {
             const bool used = dy_da[j] != 0;
             for (int k = 0; k < NA; ++k) {
-                const double a = dy_da[j] * dy_da[k];
+                const double a = dy_da[j] * dy_da[k] * w;
                 alpha[NA * j + k] += used ? a : 0.0;
             }
-            const double b = dy_sig * dy_da[j];
+            const double b = dy_sig * dy_da[j] * w;
             beta[j] += used ? b : 0.0;
         }
-        wssr += gemmi::sq(dy_sig);
+        wssr += gemmi::sq(dy_sig) * w;
     }
     // The upper half of alpha is the lower half's mirror, as in gemmi.
     for (int j = 1; j < NA; j++)
@@ -198,8 +223,8 @@ double ComputeWSSR(const FixedSolventFit &fit, std::vector<double> &kaniso) {
     for (size_t i = 0; i < fit.points.size(); ++i) {
         const auto &p = fit.points[i];
         kaniso[i] = std::exp(-0.25 * fit.b_star.r_u_r(p.hkl));
-        // FixedSolventFit::compute_value(), weight 1
-        wssr += gemmi::sq(p.get_y() - p.fcalc_abs * (float) (fit.k_overall * kaniso[i]));
+        // FixedSolventFit::compute_value()
+        wssr += gemmi::sq(p.get_y() - p.fcalc_abs * (float) (fit.k_overall * kaniso[i])) * p.get_weight();
     }
     return (double) wssr;
 }
@@ -343,10 +368,18 @@ ModelScaleReport FitModelScale(gemmi::Scaling<float> &scaling, ModelScaleBox box
     // The scale the walk stands at, which a point starts from only where the data give it no start of its own.
     double k_last = scaling.k_overall;
     gemmi::SMat33<double> b_last = scaling.b_star;
-    auto fit_point = [&scaling](FixedSolventFit &fit, PointFit &pf, double k_start,
-                                const gemmi::SMat33<double> &b_start) {
+    bool robust_pass = false;               // the fine pass of a robust fit (ModelScaleBox::robust)
+    auto fit_point = [&scaling, &robust_pass](FixedSolventFit &fit, PointFit &pf, double k_start,
+                                              const gemmi::SMat33<double> &b_start) {
         fit.Start(scaling, pf.k_sol, pf.b_sol, k_start, b_start);
         LevMarFit(fit);                      // k_overall + anisotropic B only
+        if (robust_pass) {
+            double c = 0;
+            for (int round = 0; round < ROBUST_ROUNDS; round++) {
+                RobustReweight(fit, c);
+                LevMarFit(fit);
+            }
+        }
         pf.r = RFactor(fit);
         pf.k_overall = fit.k_overall;
         pf.b_star = fit.b_star;
@@ -403,6 +436,13 @@ ModelScaleReport FitModelScale(gemmi::Scaling<float> &scaling, ModelScaleBox box
     const double k_hi2 = std::min(box.k_hi, k0 + 0.05);
     const double b_hi2 = std::min(box.b_hi, b0 + b_step);
     std::vector<PointFit> fine;
+    // A robust fit takes the fine pass robustly: each of its points, and the coarse winner again, by robust least
+    // squares - the coarse pass only finds where to look.
+    if (box.robust) {
+        robust_pass = true;
+        fitted.clear();
+        best_r = -1;
+    }
     for (double ks = std::max(box.k_lo, k0 - 0.05); ks <= k_hi2 + 1e-9; ks += 0.025)
         for (double bs = std::max(box.b_lo, b0 - b_step); bs <= b_hi2 + 1e-9; bs += b_step / 2.0)
             fine.push_back(PointFit{ks, bs});
diff --git a/image_analysis/structure_refinement/ModelScaling.h b/image_analysis/structure_refinement/ModelScaling.h
index a3d465e85..bd1434c04 100644
--- a/image_analysis/structure_refinement/ModelScaling.h
+++ b/image_analysis/structure_refinement/ModelScaling.h
@@ -31,6 +31,9 @@
 struct ModelScaleBox {
     double k_lo = 0.10, k_hi = 0.60;
     double b_lo = 10.0, b_hi = 80.0;
+    // The overall scale and anisotropic B of the points of the fine pass by robust least squares
+    // (ModelScaling.cpp, RobustReweight) rather than plain least squares on |F|: what refinement scores a model with.
+    bool robust = false;
 };
 
 // The b_sol ceiling for `model`: 80 A^2, or twice the model's mean atomic B where that is higher.
diff --git a/image_analysis/structure_refinement/Refinement.cpp b/image_analysis/structure_refinement/Refinement.cpp
index a04a80182..244fe68fe 100644
--- a/image_analysis/structure_refinement/Refinement.cpp
+++ b/image_analysis/structure_refinement/Refinement.cpp
@@ -27,6 +27,7 @@
 #include "../../common/Logger.h"
 #include "ModelFFT.h"
 #include "ModelScaling.h"
+#include "RefinementMath.h"
 #include "RefinementXray.h"
 #include "SigmaA.h"
 #ifdef JFJOCH_USE_CUDA
@@ -98,15 +99,15 @@ struct Context {
     double d_min = 0;
     RefinementReflections refl;
     std::unique_ptr<RefinementEngine> engine;
-#ifdef JFJOCH_USE_CUDA
-    std::unique_ptr<ModelScaleGPUPool> scale_gpu;
-#endif
+    std::unique_ptr<BulkSolventMask> mask;
     size_t nthreads = 1;
     int device = -1;            // the engine's GPU
     size_t engine_bytes = 0;    // and what it holds there
     XrayRowModel rows;
     std::vector<std::complex<float>> fcalc, fmask;
     double k_sol = 0, b_sol = 0;
+    double k_overall = 0;
+    gemmi::SMat33<double> b_star{0, 0, 0, 0, 0, 0};
     double b_floor = B_FLOOR;   // the cycle's lowest allowed B
     double b_ceiling = B_CEILING_MAX;
     // Riding hydrogens: per atom the parent it rides on (-1: none), and where it sits from the parent.
@@ -175,7 +176,7 @@ struct Context {
         }
         {
             const RefinementStage s2("refit.fmask");
-            engine->Fmask(Model(), fmask);
+            mask->Fmask(Model(), refl.hkl, fmask);
         }
         std::optional<RefinementStage> s3("refit.scale");
 
@@ -192,15 +193,13 @@ struct Context {
         gemmi::Scaling<float> scaling(cell, sg);
         scaling.use_solvent = true;
         scaling.prepare_points(fc, fo, &fm);
-#ifdef JFJOCH_USE_CUDA
-        if (scale_gpu)
-            scale_gpu->Fit(scaling, SolventBMax(Model()), nthreads);
-        else
-#endif
-            FitModelScale(scaling, {.b_hi = SolventBMax(Model())}, nthreads);
+        // the scale of the score: by robust least squares (ModelScaleBox::robust), on the CPU
+        FitModelScale(scaling, {.b_hi = SolventBMax(Model()), .robust = true}, nthreads);
         s3.reset();
         k_sol = scaling.k_sol;
         b_sol = scaling.b_sol;
+        k_overall = scaling.k_overall;
+        b_star = scaling.b_star;
         rows.k_total.resize(refl.Size());
         rows.fmask.resize(refl.Size());
         std::vector<std::complex<float>> fmodel(refl.Size());
@@ -217,6 +216,71 @@ struct Context {
         return ComputeRFactors(refl, fmodel);
     }
 
+    // Every term of the score of the atoms as they are, at the scale, the solvent and the likelihood parameters
+    // of the last Refit().
+    RefinementScore Score() {
+        const std::vector<std::complex<float>> fmodel = FModel();
+        RefinementScore sc;
+        sc.k_overall = k_overall;
+        sc.k_sol = k_sol;
+        sc.b_sol = b_sol;
+        sc.mask_spacing = mask->Spacing();
+        // b_star is in the reciprocal basis (exp(-h^T b_star h / 4)); in Cartesian B = O b_star O^T
+        const gemmi::Mat33 b = cell.orth.mat.multiply(b_star.as_mat33()).multiply(cell.orth.mat.transpose());
+        sc.b_iso = (b[0][0] + b[1][1] + b[2][2]) / 3.0;
+        sc.b_aniso = {b[0][0] - sc.b_iso, b[1][1] - sc.b_iso, b[2][2] - sc.b_iso, b[0][1], b[0][2], b[1][2]};
+        const RFactors r = ComputeRFactors(refl, fmodel);
+        sc.r_work = r.r_work;
+        sc.r_free = r.r_free;
+        sc.n_work = r.n_work;
+        sc.n_free = r.n_free;
+        sc.nll_work = sc.nll_free = 0;
+        double fom_sum = 0;
+        for (size_t m = 0; m < refl.Size(); m++) {
+            double dfc = 0;
+            const double fc = std::abs(fmodel[m]);
+            const double nll = refine_math::AmplitudeNegLogLikelihood<double>(refl.fo[m], fc, rows.d[m], rows.var[m],
+                                                                             refl.centric[m], dfc);
+            (refl.free[m] ? sc.nll_free : sc.nll_work) += nll;
+            // the figure of merit of the same likelihood: <cos(dphi)> = I1(X)/I0(X), X = 2 Fo D Fc / var, for an
+            // acentric reflection; tanh(X / 2) for a centric one
+            // Following Murshudov, Vagin & Dodson (1997) Acta Cryst. D53, 240-255
+            const double x = 2.0 * refl.fo[m] * rows.d[m] * fc / rows.var[m];
+            double log_i0 = 0, ratio = 0;
+            refine_math::BesselI0I1(x, log_i0, ratio);
+            fom_sum += refl.centric[m] ? std::tanh(0.5 * x) : ratio;
+        }
+        sc.fom = fom_sum / static_cast<double>(refl.Size());
+        for (int p : rider_parent)
+            sc.hydrogens += p >= 0;
+        constexpr int SHELLS = 20;
+        double s_lo = refl.inv_d2[0], s_hi = refl.inv_d2[0];
+        for (double s2 : refl.inv_d2) {
+            s_lo = std::min(s_lo, s2);
+            s_hi = std::max(s_hi, s2);
+        }
+        std::vector<std::array<double, 4>> sums(SHELLS, {0, 0, 0, 0});
+        sc.shells.resize(SHELLS);
+        for (int k = 0; k < SHELLS; k++) {
+            sc.shells[k].d_max = 1.0 / std::sqrt(s_lo + (s_hi - s_lo) * k / SHELLS);
+            sc.shells[k].d_min = 1.0 / std::sqrt(s_lo + (s_hi - s_lo) * (k + 1) / SHELLS);
+        }
+        for (size_t m = 0; m < refl.Size(); m++) {
+            const int k = std::min(SHELLS - 1, static_cast<int>((refl.inv_d2[m] - s_lo) / (s_hi - s_lo) * SHELLS));
+            const int f = refl.free[m] ? 2 : 0;
+            sums[k][f] += std::fabs(refl.fo[m] - std::abs(fmodel[m]));
+            sums[k][f + 1] += refl.fo[m];
+            (refl.free[m] ? sc.shells[k].n_free : sc.shells[k].n_work)++;
+        }
+        for (int k = 0; k < SHELLS; k++) {
+            if (sums[k][1] > 0)
+                sc.shells[k].r_work = sums[k][0] / sums[k][1];
+            if (sums[k][3] > 0)
+                sc.shells[k].r_free = sums[k][2] / sums[k][3];
+        }
+        return sc;
+    }
+
     // F_model with the cycle's scale and solvent, for the atoms as they are now.
     std::vector<std::complex<float>> FModel() {
         engine->Fcalc(Model(), fcalc);
@@ -604,7 +668,6 @@ void MakeEngine(Context &ctx, bool gpu, const RefinementSettings &settings, Logg
     const RefinementStage stage("engine");
     ctx.engine.reset();
 #ifdef JFJOCH_USE_CUDA
-    ctx.scale_gpu.reset();
     if (gpu && get_gpu_count() > 0) {
         try {
             const int device = get_gpu();
@@ -616,15 +679,12 @@ void MakeEngine(Context &ctx, bool gpu, const RefinementSettings &settings, Logg
                                                  ctx.b_ceiling, budget, &bytes, settings.nthreads);
             ctx.device = device;
             ctx.engine_bytes = bytes;
-            if (ctx.engine)
-                ctx.scale_gpu = std::make_unique<ModelScaleGPUPool>(device, ctx.refl.Size(), 1);
-            else
+            if (!ctx.engine)
                 logger.Info("Refinement: on the CPU - the GPU cannot take this case or {:.1f} GB is over its budget",
                             bytes / 1e9);
         } catch (const JFJochException &e) {
             logger.Warning("Refinement: the GPU failed ({}), on the CPU", e.what());
             ctx.engine.reset();
-            ctx.scale_gpu.reset();
         }
     }
 #else
@@ -681,6 +741,7 @@ std::unique_ptr<OwnSlot> MakeOwnSlot(const gemmi::Structure &before_topology, co
     ctx.refl = main.refl;
     ctx.nthreads = nthreads;
     ctx.b_ceiling = main.b_ceiling;
+    ctx.mask = std::make_unique<BulkSolventMask>(ctx.cell, *ctx.sg, ctx.d_min, nthreads);
     slot->topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens);
     ctx.atoms = AtomList(ctx.Model());
     ctx.rider_parent = slot->topo.rider_parent;
@@ -691,7 +752,6 @@ std::unique_ptr<OwnSlot> MakeOwnSlot(const gemmi::Structure &before_topology, co
                                                  ctx.b_ceiling, SIZE_MAX, nullptr, nthreads);
             if (!ctx.engine)
                 return nullptr;
-            ctx.scale_gpu = std::make_unique<ModelScaleGPUPool>(device, ctx.refl.Size(), 1);
         } catch (const JFJochException &) {
             return nullptr;
         }
@@ -855,6 +915,23 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
         return result;
     }
 
+    ctx.mask = std::make_unique<BulkSolventMask>(cell, sg, ctx.d_min, settings.nthreads);
+
+    // The dictionary's view of the model, and its hydrogens riding on their parents from the start, so that the
+    // model is scored as the restrained stage refines it.
+    const bool restrained = settings.restrained || settings.score_only;
+    // The weight trials' slots start from a copy of the model as it is before the topology is made.
+    const gemmi::Structure before_topology = ctx.st;
+    RefinementTopology topo;
+    if (restrained) {
+        std::optional<RefinementStage> topo_stage("topology");
+        topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens);
+        topo_stage.reset();
+        ctx.atoms = AtomList(ctx.Model());
+        ctx.rider_parent = topo.rider_parent;
+        PlaceRiders(ctx, topo);
+    }
+
     MakeEngine(ctx, settings.gpu, settings, logger);
     result.on_gpu = ctx.engine->OnGPU();
 
@@ -863,11 +940,51 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
     result.n_free = r.n_free;
     result.r_work_start = r.r_work;
     result.r_free_start = r.r_free;
-    logger.Info("Refinement: {} atoms, {} reflections to {:.2f} A ({} free), on the {}; as given R-work {:.4f}, R-free {:.4f}",
-                ctx.atoms.size(), ctx.refl.Size(), result.d_min, r.n_free, result.on_gpu ? "GPU" : "CPU", r.r_work,
-                r.r_free);
+    result.score_start = ctx.Score();
+    logger.Info("Refinement: {} atoms ({} riding hydrogens), {} reflections to {:.2f} A ({} free), on the {}; as given "
+                "R-work {:.4f}, R-free {:.4f}; k_sol {:.3f}, b_sol {:.1f}, overall B {:.1f}, mask grid {:.2f} A",
+                ctx.atoms.size(), result.score_start.hydrogens, ctx.refl.Size(), result.d_min, r.n_free,
+                result.on_gpu ? "GPU" : "CPU", r.r_work, r.r_free, ctx.k_sol, ctx.b_sol, result.score_start.b_iso,
+                result.score_start.mask_spacing);
+
+    // The maps, the blobs and the model of where the refinement ended.
+    auto finish = [&]() {
+        MakeMaps(ctx, merged, wavelength_A, result);
+        {
+            const RefinementStage stage("blobs");
+            result.blobs = FindDensityBlobs(*result.difference_map, result.difference_map_rms, ctx.st, result.d_min,
+                                            result.anomalous_map.get(), result.anomalous_map_rms);
+        }
+        result.k_sol = ctx.k_sol;
+        result.b_sol = ctx.b_sol;
+        result.model = std::make_shared<gemmi::Structure>(ctx.st);
+        result.ok = true;
+        result.seconds = std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
+    };
+    if (settings.score_only) {
+        result.score_start.hkl = ctx.refl.hkl;
+        result.score_start.fcalc = ctx.fcalc;
+        result.score_start.fmask = ctx.rows.fmask;
+        result.score_start.k_total = ctx.rows.k_total;
+        result.score_start.d = ctx.rows.d;
+        result.score_start.var = ctx.rows.var;
+        result.r_work_placement = result.r_work_start;
+        result.r_free_placement = result.r_free_start;
+        result.restrained_residues = topo.restrained_residues;
+        result.fixed_residues = topo.fixed_residues;
+        result.fixed_names = topo.fixed_names;
+        result.geometry_start = MeasureGeometry(topo.geometry, BuildRepulsions(topo, ctx.Model(), cell, sg, 0.0),
+                                                ctx.XyzB());
+        finish();
+        logger.Info("Refinement: the model scored only; geometry rmsd bonds {:.4f} A (rmsZ {:.2f}), angles {:.2f} deg "
+                    "(rmsZ {:.2f})", result.geometry_start.rmsd_bonds, result.geometry_start.rmsz_bonds,
+                    result.geometry_start.rmsd_angles, result.geometry_start.rmsz_angles);
+        return result;
+    }
 
     RefinePlacement(ctx, settings.group_b, result);
+    if (restrained)
+        PlaceRiders(ctx, topo);
     r = ctx.Refit();
     result.r_work_placement = r.r_work;
     result.r_free_placement = r.r_free;
@@ -875,18 +992,6 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
                 result.rigid_bodies, result.group_b ? "kept" : "not kept", r.r_work, r.r_free);
 
     if (settings.restrained) {
-        const gemmi::Structure before_topology = ctx.st;
-        std::optional<RefinementStage> topo_stage("topology");
-        RefinementTopology topo = BuildRefinementTopology(ctx.st, settings.riding_hydrogens);
-        topo_stage.reset();
-        // the hydrogens it added: the atoms and the engine made again for them
-        ctx.atoms = AtomList(ctx.Model());
-        ctx.rider_parent = topo.rider_parent;
-        PlaceRiders(ctx, topo);
-        if (static_cast<int>(ctx.atoms.size()) != result.atoms) {
-            MakeEngine(ctx, result.on_gpu, settings, logger);
-            result.on_gpu = ctx.engine->OnGPU();
-        }
         ctx.Refit();
         result.restrained_residues = topo.restrained_residues;
         result.fixed_residues = topo.fixed_residues;
@@ -998,17 +1103,7 @@ RefinementResult RefineModel(const gemmi::Structure &placed, const std::vector<M
         result.geometry = chosen.geometry;
     }
 
-    MakeMaps(ctx, merged, wavelength_A, result);
-    {
-        const RefinementStage stage("blobs");
-        result.blobs = FindDensityBlobs(*result.difference_map, result.difference_map_rms, ctx.st, result.d_min,
-                                        result.anomalous_map.get(), result.anomalous_map_rms);
-    }
-    result.k_sol = ctx.k_sol;
-    result.b_sol = ctx.b_sol;
-    result.model = std::make_shared<gemmi::Structure>(ctx.st);
-    result.ok = true;
-    result.seconds = std::chrono::duration<double>(std::chrono::steady_clock::now() - start).count();
+    finish();
     logger.Info("Refinement: done in {:.1f} s, {} target evaluations{}", result.seconds, ctx.evaluations,
                 result.restrained ? fmt::format("; R-work {:.4f}, R-free {:.4f}, weight {:.3g}", result.r_work,
                                                 result.r_free, result.weight)
diff --git a/image_analysis/structure_refinement/Refinement.h b/image_analysis/structure_refinement/Refinement.h
index 048439cf4..45fca1e7b 100644
--- a/image_analysis/structure_refinement/Refinement.h
+++ b/image_analysis/structure_refinement/Refinement.h
@@ -3,6 +3,9 @@
 
 #pragma once
 
+#include <array>
+#include <cmath>
+#include <complex>
 #include <memory>
 #include <string>
 #include <vector>
@@ -35,6 +38,7 @@ class Logger;
 // statistics and the map coefficients go out. Nothing it finds feeds back into processing decisions.
 struct RefinementSettings {
     bool restrained = true;          // stage 2; false stops after the placement stage
+    bool score_only = false;         // only score the model as given (RefinementResult::score_start), no stage
     bool group_b = true;             // the placement stage's B shift per residue
     int macro_cycles = 10;
     int xyz_iterations = 10;         // L-BFGS iterations per macro-cycle
@@ -49,6 +53,38 @@ struct RefinementSettings {
     size_t gpu_resident_bytes = 0;   // held on the card by what runs beside
 };
 
+// R by resolution: one of RefinementScore's shells, equal in width in 1/d^2 between the coarsest and the
+// finest reflection.
+struct RefinementShell {
+    double d_max = 0, d_min = 0;
+    int n_work = 0, n_free = 0;
+    double r_work = NAN, r_free = NAN;
+};
+
+// What scoring a model against the data gives, every term of it: the model's structure factors are
+// F_model = k_overall exp(-s^T B s / 4) (Fcalc + k_sol exp(-b_sol s^2 / 4) Fmask), the scale and the solvent
+// fitted to the working reflections; R-work and R-free; and the likelihood target, -ln L summed over the
+// working and over the free reflections, at D and the variance estimated from the free set.
+struct RefinementScore {
+    double k_overall = 0;
+    double b_iso = 0;                // the overall B: a third of the trace of B (Cartesian, A^2)
+    std::array<double, 6> b_aniso{}; // B less b_iso: B11 B22 B33 B12 B13 B23 (Cartesian, A^2)
+    double k_sol = 0, b_sol = 0;
+    double mask_spacing = 0;         // A, the grid the solvent mask is drawn on
+    double r_work = NAN, r_free = NAN;
+    int n_work = 0, n_free = 0;
+    double nll_work = NAN, nll_free = NAN;
+    double fom = NAN;                // the likelihood's mean figure of merit, all reflections
+    int hydrogens = 0;               // riding hydrogens in the model scored
+    std::vector<RefinementShell> shells;
+    // Per reflection, with settings.score_only only: the atoms' Fcalc, the solvent's k_sol exp(-b_sol s^2/4) Fmask
+    // and k_overall exp(-s^T B s / 4), in the order of `hkl`, and the likelihood's parameters.
+    std::vector<gemmi::Miller> hkl;
+    std::vector<std::complex<float>> fcalc, fmask;
+    std::vector<float> k_total;
+    std::vector<float> d, var;       // the likelihood's D and variance (epsilon beta + sigma^2)
+};
+
 // One rung of the weight ladder.
 struct RefinementWeightTrial {
     double weight = 0;               // the X-ray term's weight
@@ -66,8 +102,10 @@ struct RefinementResult {
     double d_min = 0;
     int n_work = 0, n_free = 0;
 
-    // The model as given, scaled to these reflections (hydrogens removed, B made isotropic).
+    // The model as given, scaled to these reflections (B made isotropic; the dictionary's hydrogens riding
+    // where the restrained stage runs or the model is only scored), and every term of that score.
     double r_work_start = 1, r_free_start = 1;
+    RefinementScore score_start;
     // After the placement stage.
     int rigid_bodies = 0;            // the chains moved as rigid bodies, where that was kept
     bool group_b = false;            // the B shift per residue was kept
diff --git a/image_analysis/structure_refinement/RefinementRestraints.cpp b/image_analysis/structure_refinement/RefinementRestraints.cpp
index 1f94562e3..7b3c0173c 100644
--- a/image_analysis/structure_refinement/RefinementRestraints.cpp
+++ b/image_analysis/structure_refinement/RefinementRestraints.cpp
@@ -551,17 +551,19 @@ void PlaceRidingHydrogens(RefinementTopology &topo) {
 GeometryStatistics MeasureGeometry(const RestraintSet &geometry, const RestraintSet &repulsions,
                                    const std::vector<std::array<double, 4>> &xyzb) {
     GeometryStatistics s;
-    double sb = 0, sa = 0, sc = 0;
+    double sb = 0, sa = 0, sc = 0, zb = 0, za = 0;
     int nc = 0;
     for (const RestraintRecord &r : geometry.records) {
         const RestraintKind kind = static_cast<RestraintKind>(r.kind);
         if (kind == RestraintKind::Bond) {
             const double d = RestraintValue(geometry, r, xyzb) - r.ideal;
             sb += d * d;
+            zb += d * d / (double(r.sigma) * r.sigma);
             s.bonds++;
         } else if (kind == RestraintKind::Angle) {
             const double d = RestraintValue(geometry, r, xyzb) - r.ideal;
             sa += d * d;
+            za += d * d / (double(r.sigma) * r.sigma);
             s.angles++;
         } else if (kind == RestraintKind::Chiral) {
             const double d = RestraintValue(geometry, r, xyzb) - r.ideal;
@@ -575,5 +577,7 @@ GeometryStatistics MeasureGeometry(const RestraintSet &geometry, const Restraint
     s.rmsd_bonds = s.bonds > 0 ? std::sqrt(sb / s.bonds) : 0;
     s.rmsd_angles = s.angles > 0 ? std::sqrt(sa / s.angles) : 0;
     s.rmsd_chirals = nc > 0 ? std::sqrt(sc / nc) : 0;
+    s.rmsz_bonds = s.bonds > 0 ? std::sqrt(zb / s.bonds) : 0;
+    s.rmsz_angles = s.angles > 0 ? std::sqrt(za / s.angles) : 0;
     return s;
 }
diff --git a/image_analysis/structure_refinement/RefinementRestraints.h b/image_analysis/structure_refinement/RefinementRestraints.h
index 7617e9094..969fc3a01 100644
--- a/image_analysis/structure_refinement/RefinementRestraints.h
+++ b/image_analysis/structure_refinement/RefinementRestraints.h
@@ -92,6 +92,7 @@ RestraintSet BuildRepulsions(const RefinementTopology &topo, gemmi::Model &model
 // The usual summary of a model's geometry.
 struct GeometryStatistics {
     double rmsd_bonds = 0, rmsd_angles = 0, rmsd_chirals = 0;
+    double rmsz_bonds = 0, rmsz_angles = 0;   // rms of (deviation / the dictionary's sigma)
     int bonds = 0, angles = 0;
     int clashes = 0;     // repulsions overlapping by more than 0.4 A
 };
diff --git a/image_analysis/structure_refinement/RefinementXray.cpp b/image_analysis/structure_refinement/RefinementXray.cpp
index 5ac7873f7..aa981dd0b 100644
--- a/image_analysis/structure_refinement/RefinementXray.cpp
+++ b/image_analysis/structure_refinement/RefinementXray.cpp
@@ -228,6 +228,26 @@ RefinementReflections PrepareRefinementReflections(const std::vector<MergedRefle
     return out;
 }
 
+BulkSolventMask::BulkSolventMask(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min,
+                                 size_t nthreads)
+    : nthreads_(nthreads),
+      spacing_(std::max(std::min(d_min / 3.0, MASK_SPACING_MAX), std::cbrt(cell.volume / MASK_POINTS_MAX))) {
+    grid_.unit_cell = cell;
+    grid_.spacegroup = &sg;
+    grid_.set_size_from_spacing(spacing_, gemmi::GridSizeRounding::Up);
+    orbit_leaders_ = OrbitLeaders(grid_, nthreads);
+}
+
+void BulkSolventMask::Fmask(const gemmi::Model &model, const std::vector<gemmi::Miller> &hkl,
+                            std::vector<std::complex<float>> &f) const {
+    gemmi::Grid<float> mask = grid_;
+    PutMaskOnGrid(mask, model, orbit_leaders_, nthreads_);
+    const gemmi::FPhiGrid<float> fm = MapToFPhi(mask);
+    f.clear();
+    for (const gemmi::Miller &h : hkl)
+        f.push_back(fm.get_value_by_hkl(h));
+}
+
 void EstimateLikelihoodParameters(const RefinementReflections &refl, const std::vector<std::complex<float>> &fmodel,
                                   XrayRowModel &rows, size_t nthreads) {
     const size_t n = refl.Size();
diff --git a/image_analysis/structure_refinement/RefinementXray.h b/image_analysis/structure_refinement/RefinementXray.h
index 24fd13756..336cc6d16 100644
--- a/image_analysis/structure_refinement/RefinementXray.h
+++ b/image_analysis/structure_refinement/RefinementXray.h
@@ -53,6 +53,31 @@ struct XrayRowModel {
     std::vector<float> d, var;
 };
 
+// The flat bulk-solvent model's mask (SolventMasker with the Refmac radii, PutMaskOnGrid) and its structure factors
+// at the reflections, on a grid of its own rather than the atoms' density grid: the mask's edge is set by the
+// 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. 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.
+// The flat solvent model of Fokine & Urzhumtsev (2002) Acta Cryst. D58, 1387-1392
+class BulkSolventMask {
+public:
+    static constexpr double MASK_SPACING_MAX = 0.4;   // A
+    static constexpr double MASK_POINTS_MAX = 1 << 26;
+    BulkSolventMask(const gemmi::UnitCell &cell, const gemmi::SpaceGroup &sg, double d_min, size_t nthreads);
+    // Fmask (unscaled) of `model` at `hkl`.
+    void Fmask(const gemmi::Model &model, const std::vector<gemmi::Miller> &hkl,
+               std::vector<std::complex<float>> &f) const;
+    double Spacing() const { return spacing_; }
+
+private:
+    gemmi::Grid<float> grid_;
+    std::vector<size_t> orbit_leaders_;
+    size_t nthreads_;
+    double spacing_;
+};
+
 // The ML parameters per reflection from the free reflections: D and beta fitted by maximum likelihood in
 // resolution bins of the free set and interpolated between the bins' centres; the variance of a reflection is
 // epsilon beta + sigma^2. `fmodel` is k_total (F + fmask), without D.
diff --git a/tools/rugnux_refine_test.cpp b/tools/rugnux_refine_test.cpp
index 037613063..9c367854e 100644
--- a/tools/rugnux_refine_test.cpp
+++ b/tools/rugnux_refine_test.cpp
@@ -28,6 +28,7 @@ void usage() {
                  "  -N <threads>        CPU threads (default 8)\n"
                  "  --no-gpu            the CPU engine\n"
                  "  --placement-only    rigid bodies and group B only\n"
+                 "  --score-only        score the model as given, refine nothing\n"
                  "  --no-group-b        no B shift per residue in the placement stage\n"
                  "  --no-hydrogens      no riding hydrogens\n"
                  "  --adp-ladder <a,..> the B refinement's weights, times the positions' chosen weight\n"
@@ -35,6 +36,7 @@ void usage() {
                  "  --xyz-iterations <n>, --b-iterations <n>  L-BFGS iterations per macro-cycle\n"
                  "  --ladder <a,b,...>  the weight ladder, times the gradient-balanced weight\n"
                  "  -o <file.cif>       write the refined model\n"
+                 "  --dump <file>       with --score-only: h k l |Fcalc| phi |Fmask| phi k_total D var per reflection\n"
                  "  -v                  debug lines (TIMING per refinement stage)\n";
 }
 
@@ -52,7 +54,7 @@ std::vector<double> ParseList(const std::string &s) {
 int main(int argc, char **argv) {
     RefinementSettings settings;
     settings.nthreads = 8;
-    std::string out_path;
+    std::string out_path, dump_path;
     bool verbose = false;
     std::vector<std::string> files;
     for (int i = 1; i < argc; i++) {
@@ -68,6 +70,8 @@ int main(int argc, char **argv) {
             settings.riding_hydrogens = false;
         else if (a == "--placement-only")
             settings.restrained = false;
+        else if (a == "--score-only")
+            settings.score_only = true;
         else if (a == "--cycles" && i + 1 < argc)
             settings.macro_cycles = std::stoi(argv[++i]);
         else if (a == "--xyz-iterations" && i + 1 < argc)
@@ -80,6 +84,8 @@ int main(int argc, char **argv) {
             settings.weight_ladder = ParseList(argv[++i]);
         else if (a == "-v")
             verbose = true;
+        else if (a == "--dump" && i + 1 < argc)
+            dump_path = argv[++i];
         else if (a == "-o" && i + 1 < argc)
             out_path = argv[++i];
         else if (!a.empty() && a[0] == '-') {
@@ -135,6 +141,14 @@ int main(int argc, char **argv) {
         gemmi::cif::write_cif_to_stream(f, gemmi::make_mmcif_document(*r.model));
     }
 
+    if (!dump_path.empty()) {
+        std::ofstream f(dump_path);
+        const RefinementScore &sc = r.score_start;
+        for (size_t m = 0; m < sc.hkl.size(); m++)
+            f << sc.hkl[m][0] << " " << sc.hkl[m][1] << " " << sc.hkl[m][2] << " " << std::abs(sc.fcalc[m]) << " "
+              << std::arg(sc.fcalc[m]) << " " << std::abs(sc.fmask[m]) << " " << std::arg(sc.fmask[m]) << " "
+              << sc.k_total[m] << " " << sc.d[m] << " " << sc.var[m] << "\n";
+    }
     int kinds[5] = {0, 0, 0, 0, 0};
     for (const DensityBlob &b : r.blobs)
         kinds[static_cast<int>(b.kind)]++;
@@ -158,6 +172,17 @@ int main(int argc, char **argv) {
         trials.push_back({{"weight", t.weight}, {"adp_weight", t.adp_weight}, {"r_work", t.r_work}, {"r_free", t.r_free},
                           {"rmsd_bonds", t.geometry.rmsd_bonds}, {"rmsd_angles", t.geometry.rmsd_angles}});
     j["trials"] = trials;
+    const RefinementScore &sc = r.score_start;
+    nlohmann::json shells = nlohmann::json::array();
+    for (const RefinementShell &sh : sc.shells)
+        shells.push_back({{"d_max", sh.d_max}, {"d_min", sh.d_min}, {"n_work", sh.n_work}, {"n_free", sh.n_free},
+                          {"r_work", sh.r_work}, {"r_free", sh.r_free}});
+    j["score_start"] = {{"k_overall", sc.k_overall}, {"b_iso", sc.b_iso}, {"b_aniso", sc.b_aniso}, {"k_sol", sc.k_sol},
+                        {"b_sol", sc.b_sol}, {"mask_spacing", sc.mask_spacing}, {"r_work", sc.r_work},
+                        {"r_free", sc.r_free}, {"n_work", sc.n_work}, {"n_free", sc.n_free}, {"nll_work", sc.nll_work},
+                        {"nll_free", sc.nll_free}, {"hydrogens", sc.hydrogens}, {"fom", sc.fom}, {"mean_fom", r.mean_fom},
+                        {"rmsz_bonds", r.geometry_start.rmsz_bonds}, {"rmsz_angles", r.geometry_start.rmsz_angles},
+                        {"shells", shells}};
     std::cout << "JSON " << j.dump() << std::endl;
     return EXIT_SUCCESS;
 }