Merge branch 'fix-mr-segfault' into rc175: an atom whose density at the starting radius is exactly the cutoff no longer gets a NaN radius (MR model-density overflow)

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

gemmi_gph/gemmi/dencalc.hpp

@@ -49,12 +49,15 @@ Real determine_cutoff_radius(Real x1, const ExpSum<N, Real>& precal, Real cutoff
       }
     }
   } else {
-    while (std::fabs(y2) > cutoff_level) {
+    // Local change (upstream has a while loop): a first point exactly at the cutoff, which the IT92
+    // approximation hits for some B, skipped the loop and left x1 == x2 and the radius 0/0 = NaN.
+    // Stepping out at least once gives x1 there.
+    do {
       x1 = x2;
       y1 = y2;
       x2 += 0.5f;
       y2 = precal.calculate(x2 * x2);
-    }
+    } while (std::fabs(y2) > cutoff_level);
   }
 
   return x1 + (x1 - x2) / (y1 - y2) * (cutoff_level - y1);

tests/ModelValidationTest.cpp

@@ -3,6 +3,7 @@
 
 #include <catch2/catch_all.hpp>
 
+#include <algorithm>
 #include <cmath>
 #include <cstdio>
 #include <cstring>
@@ -1003,6 +1004,42 @@ TEST_CASE("ModelValidation_ParallelGriddingMatchesGemmi", "[ModelValidation]") {
     }
 }
 
+// gemmi's radius search starts at an IT92 approximation of the radius, and for some B the density there
+// is exactly the cutoff: the search then stopped on that first point and divided 0 by 0. A nitrogen at
+// B = 277.64563 is one such B; its NaN radius gave the gridding a box it overran (the heap overflow
+// that killed molecular replacement on a large deposited model).
+TEST_CASE("ModelValidation_GriddingAnAtomWhoseRadiusStartsAtTheCutoff", "[ModelValidation]") {
+    const auto path = WriteTemp("cutoff_radius_test.pdb",
+                                ClusterPdb("CRYST1   40.000   50.000   60.000  90.00  90.00  90.00 P 1           1\n").c_str());
+    gemmi::Structure st = gemmi::read_structure_gz(path, gemmi::CoorFormat::Detect);
+    std::filesystem::remove(path);
+    for (gemmi::Chain &ch : st.models[0].chains)
+        for (gemmi::Residue &r : ch.residues)
+            for (gemmi::Atom &a : r.atoms) {
+                a.element = gemmi::Element(gemmi::El::N);
+                a.b_iso = 277.64563f;
+            }
+
+    gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
+    dc.d_min = 3.5;
+    dc.rate = 1.5;
+    dc.grid.unit_cell = st.cell;
+    dc.grid.spacegroup = st.find_spacegroup();
+    dc.blur = 0;   // the atoms' B is the B gridded
+    const auto precal = gemmi::IT92<float>::get(gemmi::El::N, 0, -1).precalculate_density_iso(277.64563f, 0.0f);
+    const float radius = dc.estimate_radius(precal, 277.64563f);
+    REQUIRE(std::isfinite(radius));
+    CHECK(radius == gemmi::it92_radius_approx(277.64563f));
+    CHECK(precal.calculate(radius * radius) == dc.cutoff);
+
+    gemmi::DensityCalculator<gemmi::IT92<float>, float> gemmi_dc = dc;
+    gemmi_dc.put_model_density_on_grid(st.models[0]);
+    PutModelDensityOnGrid(dc, st.models[0], {}, 4);
+    REQUIRE(dc.grid.data.size() == gemmi_dc.grid.data.size());
+    CHECK(std::memcmp(dc.grid.data.data(), gemmi_dc.grid.data.data(), dc.grid.data.size() * sizeof(float)) == 0);
+    CHECK(std::all_of(dc.grid.data.begin(), dc.grid.data.end(), [](float v) { return std::isfinite(v); }));
+}
+
 namespace {
     // The five groups of ModelValidation_ParallelGriddingMatchesGemmi: none, a centring, screws, a
     // cubic body centring and a cubic face centring with a diamond glide.