diff --git a/gemmi_gph/gemmi/scaling.hpp b/gemmi_gph/gemmi/scaling.hpp
index 630774c0b..519cec6e7 100644
--- a/gemmi_gph/gemmi/scaling.hpp
+++ b/gemmi_gph/gemmi/scaling.hpp
@@ -232,6 +232,11 @@ struct Scaling {
       if (p.fobs < 1 || p.fobs < p.sigma)  // skip weak reflections
         continue;
       double fcalc = std::abs(get_fcalc(p));
+      // Local change (upstream fits it): a reflection where the solvent exactly cancels the model,
+      // |Fcalc| = 0 in float, has no log ratio - its log(Fobs/0) = inf made k_overall and B NaN,
+      // and the solver that starts from them failed on a singular matrix.
+      if (fcalc == 0)
+        continue;
       double x = p.stol2;
       double y = std::log(static_cast<float>(p.fobs / fcalc));
       sx += x;
diff --git a/image_analysis/structure_refinement/ModelScaleGPU.cu b/image_analysis/structure_refinement/ModelScaleGPU.cu
index d77e4a6ed..f26a72744 100644
--- a/image_analysis/structure_refinement/ModelScaleGPU.cu
+++ b/image_analysis/structure_refinement/ModelScaleGPU.cu
@@ -71,7 +71,8 @@ __global__ void scale_sums(const int *__restrict__ hkl, const double *__restrict
         const double fo = fobs[i];
 
         if (f.mode == MODE_ISOTROPIC) {
-            if (!strong[i])
+            // fit_isotropic_b_approximately() (local change in gemmi_gph): no log ratio at |F| = 0
+            if (!strong[i] || f_abs == 0)
                 continue;
             const double x = stol2[i];
             const double y = logf(static_cast<float>(fo / f_abs));
diff --git a/image_analysis/structure_refinement/ModelScaling.cpp b/image_analysis/structure_refinement/ModelScaling.cpp
index d7b941c5a..6a099189b 100644
--- a/image_analysis/structure_refinement/ModelScaling.cpp
+++ b/image_analysis/structure_refinement/ModelScaling.cpp
@@ -65,6 +65,8 @@ struct FixedSolventFit {
             if (p.fobs < 1 || p.fobs < p.sigma)  // skip weak reflections
                 continue;
             double fcalc_abs = points[i].fcalc_abs;
+            if (fcalc_abs == 0)  // local change in gemmi_gph: no log ratio
+                continue;
             double x = p.stol2;
             double y = std::log(static_cast<float>(p.fobs / fcalc_abs));
             sx += x;
@@ -228,6 +230,11 @@ void LevMarFit(FixedSolventFit &target) {
     for (;;) {
         if (eval_limit > 0 && eval_count >= eval_limit)
             break;
+        // Not in gemmi: a fit that has walked to where k_overall exp(-B s^2) |F| underflows at every
+        // reflection has no step left to take - the squares in alpha are all 0 while beta keeps a
+        // residue near 1e-210, and jordan_solve() throws on that. It ends here, at its best so far.
+        if (alpha[0] == 0)
+            break;
 
         // prepare next parameters -> temp_beta
         std::copy(alpha, alpha + NA * NA, temp_alpha);
diff --git a/tests/ModelScaleGPUTest.cpp b/tests/ModelScaleGPUTest.cpp
index 86c167f83..67a52067f 100644
--- a/tests/ModelScaleGPUTest.cpp
+++ b/tests/ModelScaleGPUTest.cpp
@@ -382,4 +382,24 @@ TEST_CASE("ModelScaleGPUPool_MatchesFitModelScale", "[ModelValidation][gpu]") {
     }
 }
 
+// ModelScaling fits past a reflection the solvent cancels exactly (ModelScalingTest.cpp), on the card.
+TEST_CASE("ModelScaleGPUPool_FitsPastAReflectionTheSolventCancelsExactly", "[ModelValidation][gpu]") {
+    if (get_gpu_count() == 0)
+        SKIP("No GPU");
+    gemmi::UnitCell cell(72, 44, 51, 90, 112, 90);
+    gemmi::Scaling<float> scaling = MakeScaling(cell, *gemmi::find_spacegroup_by_name("C 1 2 1"), 2.5);
+    // exp(-b_sol s^2) = 1 at s = 0, so at the grid's k_sol = 0.10 the sum is -0.1f + 0.1f = 0 for every b_sol
+    auto &p = *std::find_if(scaling.points.begin(), scaling.points.end(),
+                            [](const auto &q) { return q.fobs >= 1 && q.fobs >= q.sigma; });
+    p.stol2 = 0;
+    p.fcmol = {-0.1f, 0.f};
+    p.fmask = {1.f, 0.f};
+
+    ModelScaleGPUPool pool(0, scaling.points.size(), 1);
+    ModelScaleReport report;
+    REQUIRE_NOTHROW(report = pool.Fit(scaling, 80.0, 2));
+    CHECK(std::isfinite(report.r_work_fit));
+    CHECK(std::isfinite(scaling.k_overall));
+}
+
 #endif
diff --git a/tests/ModelScalingTest.cpp b/tests/ModelScalingTest.cpp
index b5b5247f7..3898b033d 100644
--- a/tests/ModelScalingTest.cpp
+++ b/tests/ModelScalingTest.cpp
@@ -271,3 +271,36 @@ TEST_CASE("ModelScaling fits a grid point exactly as gemmi's Scaling does") {
         CHECK(ours.calculate_r_factor() == gemmis.calculate_r_factor());
     }
 }
+
+// A reflection where the solvent term cancels the model exactly - |Fcalc + k_sol exp(-b_sol s^2) Fmask| = 0
+// in float - at some grid points. Its log(Fobs/|Fcalc|) was inf, which made the starting scale there NaN and
+// the solver fail on a singular matrix. It has no log ratio, so the starting fit leaves it out.
+TEST_CASE("ModelScaling fits past a reflection the solvent cancels exactly") {
+    const gemmi::UnitCell cell = Cell();
+    const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 21 21 21");
+    REQUIRE(sg != nullptr);
+    gemmi::Scaling<float> scaling(cell, sg);
+    scaling.use_solvent = true;
+    MakePoints(scaling, cell, /*k_overall=*/0.5, /*b_iso=*/20.0,
+               /*k_sol=*/0.35, /*b_sol=*/45.0, /*d_min=*/2.0, /*d_max=*/50.0);
+    // exp(-b_sol s^2) = 1 at s = 0, so at the grid's k_sol = 0.10 the sum is -0.1f + 0.1f = 0 for every b_sol
+    auto &p = scaling.points[0];
+    p.stol2 = 0;
+    p.fcmol = {-0.1f, 0.f};
+    p.fmask = {1.f, 0.f};
+    REQUIRE(p.fobs >= 1);
+
+    gemmi::Scaling<float> gemmis = scaling;
+    ModelScaleReport report;
+    REQUIRE_NOTHROW(report = FitModelScale(scaling));
+    CHECK(std::isfinite(report.r_work_fit));
+
+    gemmis.fix_k_sol = true;
+    gemmis.fix_b_sol = true;
+    gemmis.k_sol = 0.10;
+    gemmis.b_sol = 10.0;
+    REQUIRE(std::abs(gemmis.get_fcalc(gemmis.points[0])) == 0);
+    gemmis.fit_isotropic_b_approximately();
+    CHECK(std::isfinite(gemmis.k_overall));
+    REQUIRE_NOTHROW(gemmis.fit_parameters());
+}