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

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.