diff --git a/image_analysis/structure_refinement/ModelScaling.cpp b/image_analysis/structure_refinement/ModelScaling.cpp
index 8ce43a4e2..ba9c42df6 100644
--- a/image_analysis/structure_refinement/ModelScaling.cpp
+++ b/image_analysis/structure_refinement/ModelScaling.cpp
@@ -105,11 +105,11 @@ struct FixedSolventFit {
     double compute_value(const Point &p) const {
         return p.fcalc_abs * (float) get_overall_scale_factor(p.hkl);
     }
-    // For NA parameters: the B components are NA - 1 constraint rows.
+    // For NA parameters: the B components are NA - 1 constraint rows. `kaniso` is
+    // exp(-0.25 * b_star.r_u_r(h)) at the current parameters, as compute_value() takes it too.
     template <int NA>
-    double compute_value_and_derivatives(const Point &p, double *dy_da) const {
+    double compute_value_and_derivatives(const Point &p, double kaniso, double *dy_da) const {
         gemmi::Vec3 h(p.hkl);
-        double kaniso = std::exp(-0.25 * b_star.r_u_r(h));
         double fcalc_abs = p.fcalc_abs;
         int n = 1;
         double fe = fcalc_abs * kaniso;
@@ -149,13 +149,18 @@ double RFactor(const FixedSolventFit &fit) {
 // adding +0 would change, since it starts at +0 and a sum is only -0 when both its terms are. And
 // where gemmi fills the lower half of alpha, the whole square is summed here, a rectangle the compiler
 // can vectorise; the lower half is gemmi's, and the upper is replaced by it at the end, as in gemmi.
+//
+// `kaniso`: each reflection's exp(-0.25 * b_star.r_u_r(h)) at the current parameters where ComputeWSSR()
+// has just taken them there (nullptr: computed here) - the same value, which need not be taken twice.
 template <int NA>
-double ComputeLMMatrices(const FixedSolventFit &fit, double *alpha_out, double *beta_out) {
+double ComputeLMMatrices(const FixedSolventFit &fit, const double *kaniso, double *alpha_out, double *beta_out) {
     long double wssr = 0; // long double here notably increases the accuracy
     double alpha[NA * NA] = {}, beta[NA] = {};
     double dy_da[NA];
-    for (const auto &p : fit.points) {
-        double y = fit.compute_value_and_derivatives<NA>(p, dy_da);
+    for (size_t i = 0; i < fit.points.size(); ++i) {
+        const auto &p = fit.points[i];
+        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;
         for (int j = 0; j != NA; ++j) {
             const bool used = dy_da[j] != 0;
@@ -177,8 +182,21 @@ double ComputeLMMatrices(const FixedSolventFit &fit, double *alpha_out, double *
     return (double) wssr;
 }
 
-// gemmi's LevMar::fit() (levmar.hpp) with its default settings, on ComputeLMMatrices<NA>() - the same
-// iterations to the same answer, bit for bit.
+// gemmi's compute_wssr() (levmar.hpp) on a FixedSolventFit, keeping each reflection's
+// exp(-0.25 * b_star.r_u_r(h)) in `kaniso` for ComputeLMMatrices() at the same parameters.
+double ComputeWSSR(const FixedSolventFit &fit, std::vector<double> &kaniso) {
+    long double wssr = 0; // long double here notably increases the accuracy
+    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]));
+    }
+    return (double) wssr;
+}
+
+// gemmi's LevMar::fit() (levmar.hpp) with its default settings, on ComputeLMMatrices<NA>() and
+// ComputeWSSR() - the same iterations to the same answer, bit for bit.
 template <int NA>
 void LevMarFit(FixedSolventFit &target) {
     const int eval_limit = 100;
@@ -194,7 +212,8 @@ void LevMarFit(FixedSolventFit &target) {
     double alpha[NA * NA], beta[NA], temp_alpha[NA * NA];
     std::vector<double> temp_beta(NA);
 
-    const double initial_wssr = ComputeLMMatrices<NA>(target, alpha, beta);
+    std::vector<double> kaniso(target.points.size());
+    const double initial_wssr = ComputeLMMatrices<NA>(target, nullptr, alpha, beta);
     double wssr = initial_wssr;
 
     int small_change_counter = 0;
@@ -219,7 +238,7 @@ void LevMarFit(FixedSolventFit &target) {
             temp_beta[i] += best_a[i];
 
         target.set_parameters(temp_beta);
-        double new_wssr = gemmi::compute_wssr(target);
+        double new_wssr = ComputeWSSR(target, kaniso);
         ++eval_count;
         if (new_wssr < wssr) {
             double rel_change = (wssr - new_wssr) / wssr;
@@ -235,7 +254,7 @@ void LevMarFit(FixedSolventFit &target) {
             } else {
                 small_change_counter = 0;
             }
-            ComputeLMMatrices<NA>(target, alpha, beta);
+            ComputeLMMatrices<NA>(target, kaniso.data(), alpha, beta);   // at the parameters just tried
             ++eval_count;
             lambda *= lambda_down_factor;
         } else { // worse fitting