RESULT-CHANGING (ulp level, CUDA builds): the CUDA model density is the CPU walk's, bit for bit

ModelDensityGPU's gather had arithmetic of its own: per-thread offsets in Cartesian, the atom's
nearest image to the brick, r2 = x x + y y + z z, the radius squared in float. It now does what
UseDensityBoxesByPlane + AtomDensity (ModelGrid.cpp) and ModelDensityMetal.metal do, in their order:
- each brick's atoms carry their box (c - d, wrapped start, 2 d); a point is in the box by a
  subtraction, and f = centre - (u, v, w) in the box's own unwrapped index;
- x = m00 fx + (m01 fy + m02 fz), y = m11 fy + m12 fz, z = m22 fz with gemmi's orth_n
  (orth * (1/n), as Grid sets it) in float; r2 = (y y + z z) + x x against max_r2 = float(r^2)
  rounded from the double radius (new ModelDensityAtom field);
- gemmi's expapprox, the five Gaussians in order, atoms in model order per point. Device code is
  already compiled with --fmad=false (top-level CMakeLists), so nothing is contracted.
- The centres are no longer wrapped into the cell (ModelStructureFactorsGPU, RigidBodyGPU's
  place_kernel): float(fractional * n), as MakeDensityBox holds them; the mask keeps wrapped
  fractional coordinates. The 'narrow' nearest-image branch is gone: Supports() already leaves
  only boxes narrower than the cell.
The target is PutModelDensityOnGrid() with no orbits to symmetrize (one copy of the content, which
is what the gather grids; symmetry stays composed in reciprocal space). A rounding tie in the box
centre c (float on the GPU, double on the CPU) changes no point: either box covers the sphere with
half a spacing to spare. New ModelDensityGPU_BitIdenticalToCPU ([ModelValidation][gpu]) memcmps the
grids: triclinic P1 and P 21 21 21 cells, 3000 atoms over -0.5..1.5 of the cell, a quarter
anisotropic - 0 points differ.

RTX 5080, alternating base (8d199cceb) / this, under the GPU lock:
- gather kernel (nsys): 4nwv --model median 2.06 -> 2.55 ms, total 0.79 -> 0.97 s;
  9s02 --mr median 0.30 -> 0.37 ms, total 0.94 -> 1.17 s (+~25%, the box test and per-point
  index arithmetic).
- wall, two each: 4nwv --model 76.0/68.9 -> 74.4/71.0 s, --mr 71.3/71.1 -> 70.9/71.6 s;
  9s02 --model 25.1/22.6 -> 23.3/23.5 s, --mr 23.3/23.6 -> 23.3/23.2 s; 5reo --model
  6.2/5.5 -> 5.5/5.7 s. Model validation 4nwv 12.4/12.3 -> 12.3/11.7 s; MR 4nwv 15.4/15.2 ->
  15.4/15.1 s, 9s02 2.3/2.4 -> 2.4/2.4 s: no change beyond noise.
- outputs (each build run-to-run identical): p.hkl identical on all; R-work/R-free and every
  decision unchanged (4nwv validation 0.2982/0.2969, 9s02 MR AGREES P 21 21 2 0.2503/0.2517,
  4nwv MR UNDECIDED I 2 2 2 / I 21 21 21). Moved: MODEL_FIT_SIGMA 9s02 99.21 -> 101.68, 5reo
  103.96 -> 104.98; MR placement Z 9s02 142.1 -> 145.6, 4nwv 46.7 -> 46.2; 4nwv anomalous site
  4.01 -> 4.02 sigma; 4nwv placed model 1322/125462 lines in the last digit (2fo-fc max |d| 0.04
  rms, CC 0.99999994); 9s02/5reo models identical, maps to 2e-5 rms.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01PXc1zD1yxHqaR1DGdCTHL4

image_analysis/structure_refinement/ModelStructureFactorsGPU.cpp

@@ -77,7 +77,9 @@ void ModelDensityAtoms(const gemmi::Model &model, const gemmi::DensityCalculator
                 if (!a.aniso) {
                     const CReal b = static_cast<CReal>(atom.b_iso + dc.blur);
                     const auto precal = coef.precalculate_density_iso(b, addend);
-                    a.radius = dc.estimate_radius(precal, b);
+                    const double radius = dc.estimate_radius(precal, b);
+                    a.radius = static_cast<float>(radius);
+                    a.max_r2 = static_cast<float>(radius * radius);
                     for (int k = 0; k < 5; k++) {
                         a.a[k] = precal.a[k];
                         a.b[k][0] = precal.b[k];
@@ -85,7 +87,9 @@ void ModelDensityAtoms(const gemmi::Model &model, const gemmi::DensityCalculator
                 } else {
                     const auto aniso_b = atom.aniso.scaled(CReal(gemmi::u_to_b())).added_kI(CReal(dc.blur));
                     const CReal b_max = std::max(std::max(aniso_b.u11, aniso_b.u22), aniso_b.u33);
-                    a.radius = static_cast<float>(dc.estimate_radius(coef.precalculate_density_iso(b_max, addend), b_max));
+                    const double radius = dc.estimate_radius(coef.precalculate_density_iso(b_max, addend), b_max);
+                    a.radius = static_cast<float>(radius);
+                    a.max_r2 = static_cast<float>(radius * radius);
                     const auto precal = coef.precalculate_density_aniso_b(aniso_b, addend);
                     for (int k = 0; k < 5; k++) {
                         a.a[k] = precal.a[k];
@@ -191,7 +195,8 @@ void ModelStructureFactorsGPU::Compute(const gemmi::Model &model, gemmi::AsuData
     std::vector<float> mask_radius;
     ModelDensityAtoms(model, dc_, atoms, mask_atom, mask_radius);
 
-    // Each atom in grid units and, for the mask, fractional - both wrapped into the cell.
+    // Each atom in grid units, as the CPU's density walk centres it (MakeDensityBox, ModelGrid.cpp), and, for
+    // the mask, fractional and wrapped into the cell.
     const int n3[3] = {setup_.grid.nu, setup_.grid.nv, setup_.grid.nw};
     std::vector<std::array<float, 4>> pos;
     std::vector<std::array<double, 3>> frac;
@@ -202,12 +207,10 @@ void ModelStructureFactorsGPU::Compute(const gemmi::Model &model, gemmi::AsuData
                 std::array<double, 3> f{f0.x, f0.y, f0.z};
                 std::array<float, 4> g{};
                 for (int k = 0; k < 3; k++) {
+                    g[k] = static_cast<float>(f[k] * n3[k]);
                     f[k] -= std::floor(f[k]);
                     if (f[k] >= 1.0)   // a tiny negative coordinate wraps to 1 - epsilon, which rounds to 1
                         f[k] = 0.0;
-                    g[k] = static_cast<float>(f[k] * n3[k]);
-                    if (g[k] >= n3[k])
-                        g[k] -= n3[k];
                 }
                 frac.push_back(f);
                 pos.push_back(g);