From 29e4887efa0b8582473a96fa0a5c6259ce6044a7 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Tue, 13 Oct 2026 14:28:51 +0200
Subject: [PATCH] gemmi dencalc: an atom whose density at the starting radius
 is exactly the cutoff gets a radius, not NaN

determine_cutoff_radius() starts at the IT92 approximation of the radius. When the density there is
exactly the float cutoff (1e-5f) - which happens for a few B per element, e.g. a nitrogen at
B + blur = 277.64563 - neither loop runs, x1 == x2, and the interpolation returns 0/0 = NaN. Upstream
gemmi master has the same code.

A NaN radius cast to int gave MakeAtomBox() INT_MIN half-widths, so UseAtomBoxesByPlane() sized its
point buffer from a wrapped product (one point) and then wrote past it: the heap overflow (SIGSEGV, or
a crash later in FFTW) that killed molecular replacement after the rigid body of a large deposited
model, turned to imitate another crystal form, had failed (2 of 73476 atoms, in the CPU rigid body of
an MR placement). The GPU density path took the same NaN radius into its box sizes.

The outward search now steps at least once, so that case returns the point at the cutoff. Every other
B gives the same radius bit for bit (checked over ten elements, every 8th float B in 2-600 A^2: 3 NaN
before, none after, no finite radius changed).

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01EizBKhTqrAYago9KJG3dA3
---
 gemmi_gph/gemmi/dencalc.hpp   |  7 +++++--
 tests/ModelValidationTest.cpp | 37 +++++++++++++++++++++++++++++++++++
 2 files changed, 42 insertions(+), 2 deletions(-)

diff --git a/gemmi_gph/gemmi/dencalc.hpp b/gemmi_gph/gemmi/dencalc.hpp
index f5d9217ba..184f44d3b 100644
--- a/gemmi_gph/gemmi/dencalc.hpp
+++ b/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);
diff --git a/tests/ModelValidationTest.cpp b/tests/ModelValidationTest.cpp
index 5bac11561..2af466c49 100644
--- a/tests/ModelValidationTest.cpp
+++ b/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.