Merge cuda-density-align into rc176: the CUDA model density is the CPU walk's bit for bit (RESULT-CHANGING at the last bit on CUDA builds)

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

image_analysis/structure_refinement/ModelDensityGPU.cu

@@ -112,31 +112,41 @@ __global__ void brick_ranges_kernel(const int *key, int npairs, int *start, int
         end[k] = i + 1;
 }
 
-// The density of the model on the grid, as gemmi's do_add_atom_density_to_grid() puts it there: every
-// point adds, in model order, each atom whose sphere it is inside. One block per brick; the brick's atoms
-// are staged in shared memory with the position of their image nearest to the brick, relative to the
-// brick's first point. Where the cell is wider than an atom's box plus a brick, that is the only image
-// that reaches any point of the brick; on a narrower cell each point finds its own nearest image.
-__global__ void gather_kernel(const ModelDensityAtom *atoms, const float4 *pos, const int *pair_atom,
+// An atom's box as the gather tests it: its first point, unwrapped (c - d) and wrapped into the grid, and its
+// width less one (2 d), along each axis; and its centre in grid units.
+struct GatherBox {
+    int lo[3], first[3], span[3];
+    float nctr[3];
+};
+
+// The point p's offset in the box along an axis of n points, 0 to span where it is in it.
+__device__ __forceinline__ int box_offset(int p, int first, int n) {
+    const int m = p - first;
+    return m < 0 ? m + n : m;
+}
+
+// The density of the model on the grid, as UseDensityBoxesByPlane() and AtomDensity (ModelGrid.cpp) put
+// it there, bit for bit: every point adds, in model order, each atom whose box holds it and whose radius
+// reaches it, with the CPU's single-precision arithmetic in the CPU's order (no a*b+c is contracted:
+// CMakeLists.txt compiles device code with --fmad=false). One block per brick; the brick's atoms are
+// staged in shared memory. Supports() leaves only boxes narrower than the cell, so a point is in an atom's
+// box at most once, at one index of the box's own unwrapped range.
+__global__ void gather_kernel(const ModelDensityAtom *atoms, const float4 *pos, const int *box, const int *pair_atom,
                               const int *brick_start, const int *brick_end, ModelDensityGPU::Geometry g,
                               float *grid) {
     const int nbu = (g.nu + BRICK - 1) / BRICK, nbv = (g.nv + BRICK - 1) / BRICK;
     const int brick = blockIdx.x;
     const int s = brick_start[brick], e = brick_end[brick];   // s == e: an empty brick, written as zeros
-    const int u0 = (brick % nbu) * BRICK, v0 = (brick / nbu % nbv) * BRICK, w0 = (brick / (nbu * nbv)) * BRICK;
-    const int u = u0 + threadIdx.x, v = v0 + threadIdx.y, w = w0 + threadIdx.z;
+    const int u = (brick % nbu) * BRICK + threadIdx.x, v = (brick / nbu % nbv) * BRICK + threadIdx.y,
+              w = (brick / (nbu * nbv)) * BRICK + threadIdx.z;
     const bool inside = u < g.nu && v < g.nv && w < g.nw;
     const int tid = threadIdx.x + BRICK * (threadIdx.y + BRICK * threadIdx.z);
-    // This point relative to the brick's first one, in Cartesian coordinates.
-    const float tx = g.orth_n[0] * threadIdx.x + g.orth_n[1] * threadIdx.y + g.orth_n[2] * threadIdx.z;
-    const float ty = g.orth_n[3] * threadIdx.x + g.orth_n[4] * threadIdx.y + g.orth_n[5] * threadIdx.z;
-    const float tz = g.orth_n[6] * threadIdx.x + g.orth_n[7] * threadIdx.y + g.orth_n[8] * threadIdx.z;
+    const int n3[3] = {g.nu, g.nv, g.nw};
 
     constexpr int WORDS = sizeof(ModelDensityAtom) / sizeof(int);
     __shared__ int sh_index[GATHER_TILE];
     __shared__ ModelDensityAtom sh_atom[GATHER_TILE];
-    __shared__ float3 sh_centre[GATHER_TILE];   // Cartesian
-    __shared__ float3 sh_grid[GATHER_TILE];     // the same, in grid units
+    __shared__ GatherBox sh_box[GATHER_TILE];
     float acc = 0.0f;
     for (int c = s; c < e; c += GATHER_TILE) {
         const int m = min(GATHER_TILE, e - c);
@@ -144,16 +154,16 @@ __global__ void gather_kernel(const ModelDensityAtom *atoms, const float4 *pos,
         if (tid < m) {
             const int a = pair_atom[c + tid];
             sh_index[tid] = a;
-            // The atom's image nearest to the brick's centre, in grid units from its first point.
             const float4 p = pos[a];
-            float fx = p.x - u0, fy = p.y - v0, fz = p.z - w0;
-            fx -= g.nu * rintf((fx - 0.5f * (BRICK - 1)) / g.nu);
-            fy -= g.nv * rintf((fy - 0.5f * (BRICK - 1)) / g.nv);
-            fz -= g.nw * rintf((fz - 0.5f * (BRICK - 1)) / g.nw);
-            sh_grid[tid] = make_float3(fx, fy, fz);
-            sh_centre[tid] = make_float3(g.orth_n[0] * fx + g.orth_n[1] * fy + g.orth_n[2] * fz,
-                                         g.orth_n[3] * fx + g.orth_n[4] * fy + g.orth_n[5] * fz,
-                                         g.orth_n[6] * fx + g.orth_n[7] * fy + g.orth_n[8] * fz);
+            const float nctr[3] = {p.x, p.y, p.z};
+            GatherBox &b = sh_box[tid];
+            for (int k = 0; k < 3; k++) {
+                const int d = box[3 * a + k];
+                b.lo[k] = __float2int_rn(nctr[k]) - d;
+                b.first[k] = imod(b.lo[k], n3[k]);
+                b.span[k] = 2 * d;
+                b.nctr[k] = nctr[k];
+            }
         }
         __syncthreads();
         for (int i = tid; i < m * WORDS; i += BRICK * BRICK * BRICK)
@@ -162,19 +172,22 @@ __global__ void gather_kernel(const ModelDensityAtom *atoms, const float4 *pos,
         if (!inside)
             continue;
         for (int k = 0; k < m; k++) {
+            const GatherBox &b = sh_box[k];
+            const int mu = box_offset(u, b.first[0], g.nu), mv = box_offset(v, b.first[1], g.nv),
+                      mw = box_offset(w, b.first[2], g.nw);
+            if (mu > b.span[0] || mv > b.span[1] || mw > b.span[2])
+                continue;
             const ModelDensityAtom &at = sh_atom[k];
-            float x = sh_centre[k].x - tx, y = sh_centre[k].y - ty, z = sh_centre[k].z - tz;
-            if (g.narrow) {
-                float fx = sh_grid[k].x - threadIdx.x, fy = sh_grid[k].y - threadIdx.y, fz = sh_grid[k].z - threadIdx.z;
-                fx -= g.nu * rintf(fx / g.nu);
-                fy -= g.nv * rintf(fy / g.nv);
-                fz -= g.nw * rintf(fz / g.nw);
-                x = g.orth_n[0] * fx + g.orth_n[1] * fy + g.orth_n[2] * fz;
-                y = g.orth_n[3] * fx + g.orth_n[4] * fy + g.orth_n[5] * fz;
-                z = g.orth_n[6] * fx + g.orth_n[7] * fy + g.orth_n[8] * fz;
-            }
-            const float r2 = x * x + y * y + z * z;
-            if (r2 > at.radius * at.radius)
+            // f = centre - (u, v, w) in the box's own (unwrapped) indices, then in A, as the CPU walks it
+            const float fx = b.nctr[0] - static_cast<float>(b.lo[0] + mu);
+            const float fy = b.nctr[1] - static_cast<float>(b.lo[1] + mv);
+            const float fz = b.nctr[2] - static_cast<float>(b.lo[2] + mw);
+            const float z = g.m22 * fz;
+            const float y = g.m11 * fy + g.m12 * fz;
+            const float yz = y * y + z * z;
+            const float x = g.m00 * fx + (g.m01 * fy + g.m02 * fz);
+            const float r2 = yz + x * x;
+            if (r2 > at.max_r2)
                 continue;
             float density = 0.0f;
             if (!at.aniso) {
@@ -182,8 +195,10 @@ __global__ void gather_kernel(const ModelDensityAtom *atoms, const float4 *pos,
                     density += at.a[q] * expapprox(fmaxf(at.b[q][0] * r2, -88.f));
             } else {
                 for (int q = 0; q < 5; q++) {
-                    const float *b = at.b[q];
-                    const float rur = x * x * b[0] + y * y * b[1] + z * z * b[2] + 2 * (x * y * b[3] + x * z * b[4] + y * z * b[5]);
+                    const float *bq = at.b[q];
+                    // gemmi's SMat33::r_u_r(), in its order
+                    const float rur = x * x * bq[0] + y * y * bq[1] + z * z * bq[2] +
+                                      2 * (x * y * bq[3] + x * z * bq[4] + y * z * bq[5]);
                     density += at.a[q] * expapprox(fmaxf(rur, -88.f));
                 }
             }
@@ -288,20 +303,21 @@ void ModelDensityGPU::SetAtoms(const ModelDensityGrid &grid, const std::vector<M
     if (atom_list.size() > max_atoms || Bricks(grid.nu, grid.nv, grid.nw) > max_bricks)
         throw JFJochException(JFJochExceptionCategory::GPUCUDAError, "model density: grid or atoms over the reserve");
     // An atom whose box is wider than the cell would reach a point through more than one image, which
-    // gemmi's box walk adds and the gather's nearest image does not; such a cell is left to the CPU.
+    // gemmi's box walk adds and the gather's box test does not; such a cell is left to the CPU.
     if (!Supports(grid, atom_list))
         throw JFJochException(JFJochExceptionCategory::GPUCUDAError, "model density: an atom wider than the cell");
     geom.nu = grid.nu;
     geom.nv = grid.nv;
     geom.nw = grid.nw;
-    for (int i = 0; i < 3; i++)
-        for (int j = 0; j < 3; j++)
-            geom.orth_n[3 * i + j] = static_cast<float>(grid.orth[3 * i + j] / (j == 0 ? grid.nu : j == 1 ? grid.nv : grid.nw));
+    // As gemmi's Grid sets orth_n (multiply_by_diagonal with 1/n) and ModelDensityGeometry() rounds it.
+    const double inv_n[3] = {1.0 / grid.nu, 1.0 / grid.nv, 1.0 / grid.nw};
+    geom.m00 = static_cast<float>(grid.orth[0] * inv_n[0]);
+    geom.m01 = static_cast<float>(grid.orth[1] * inv_n[1]);
+    geom.m02 = static_cast<float>(grid.orth[2] * inv_n[2]);
+    geom.m11 = static_cast<float>(grid.orth[4] * inv_n[1]);
+    geom.m12 = static_cast<float>(grid.orth[5] * inv_n[2]);
+    geom.m22 = static_cast<float>(grid.orth[8] * inv_n[2]);
     const std::vector<int> b = Boxes(grid, atom_list);
-    geom.narrow = 0;
-    for (size_t i = 0; i < b.size(); i++)
-        if (2 * b[i] + BRICK > (i % 3 == 0 ? grid.nu : i % 3 == 1 ? grid.nv : grid.nw))
-            geom.narrow = 1;
     n_atoms = static_cast<int>(atom_list.size());
     if (n_atoms > 0) {
         cuda_err(cudaMemcpyAsync(atoms, atom_list.data(), atom_list.size() * sizeof(ModelDensityAtom),
@@ -344,7 +360,7 @@ void ModelDensityGPU::Compute(const float4 *d_pos, float *d_grid) {
         brick_ranges_kernel<<<blocks(npairs), THREADS, 0, stream>>>(key_sorted, npairs, brick_start, brick_end);
         cuda_err(cudaGetLastError());
     }
-    gather_kernel<<<nb, dim3(BRICK, BRICK, BRICK), 0, stream>>>(atoms, d_pos, value_sorted, brick_start, brick_end,
+    gather_kernel<<<nb, dim3(BRICK, BRICK, BRICK), 0, stream>>>(atoms, d_pos, box, value_sorted, brick_start, brick_end,
                                                                 geom, d_grid);
     cuda_err(cudaGetLastError());
 }

image_analysis/structure_refinement/ModelDensityGPU.h

@@ -9,8 +9,11 @@
 // with it, and compose the crystal's symmetry in reciprocal space (SymmetryComposition, RigidBodyRefine.h).
 //
 // A gather, not a scatter: the grid is cut into bricks of BRICK^3 points, each brick is one block, and
-// every point adds, in model order, each atom whose sphere it is inside. Nothing is summed with a
-// floating-point atomic, so the grid is the same, bit for bit, on every run.
+// every point adds, in model order, each atom whose box holds it and whose radius reaches it. Nothing is
+// summed with a floating-point atomic. The arithmetic is the CPU walk's (UseDensityBoxesByPlane and
+// AtomDensity, ModelGrid.cpp; ModelDensityMetal.metal spells out the same), so the grid is
+// PutModelDensityOnGrid()'s before its symmetrize, bit for bit, for the same atoms and centres -
+// tests/ModelDensityGPUTest.cpp holds it to that.
 
 #include <cstddef>
 #include <vector>
@@ -26,6 +29,7 @@ struct ModelDensityAtom {
     float b[5][6];      // isotropic: b[k][0] multiplies r^2; anisotropic: the matrix, u11 u22 u33 u12 u13 u23
     float occ;
     float radius;
+    float max_r2;       // the radius squared, as PutModelDensityOnGrid() rounds it: float(radius^2) in double
     int aniso;
 };
 
@@ -56,16 +60,17 @@ public:
     // The grid and the atoms' densities. Throws if they are over the capacity or if !Supports().
     void SetAtoms(const ModelDensityGrid &grid, const std::vector<ModelDensityAtom> &atoms);
 
-    // The density into d_grid (nu * nv * nw floats, every point written). d_pos: each atom's position in
-    // grid units - fractional times the grid size, wrapped into [0, n) - in the order SetAtoms() got them;
-    // w is not read. Waits on the stream once, for the number of (brick, atom) pairs.
+    // The density into d_grid (nu * nv * nw floats, every point written). d_pos: each atom's centre in
+    // grid units - its fractional position times the grid size, as the CPU's DensityBox holds it, not
+    // wrapped into the cell - in the order SetAtoms() got them; w is not read. Waits on the stream once,
+    // for the number of (brick, atom) pairs.
     void Compute(const float4 *d_pos, float *d_grid);
 
     // What the gather kernel needs of the grid.
     struct Geometry {
         int nu, nv, nw;
-        float orth_n[9];        // orth * diag(1/nu, 1/nv, 1/nw), row-major: grid offset -> Cartesian
-        int narrow;             // the cell is too small for one image of an atom to serve a whole brick
+        // gemmi's orth_n, orth * diag(1/nu, 1/nv, 1/nw), upper triangular, in single precision
+        float m00, m01, m02, m11, m12, m22;
     };
 
 private:

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);

image_analysis/structure_refinement/RigidBodyGPU.cu

@@ -85,8 +85,9 @@ __device__ __forceinline__ double r_u_r(const double *u, const int *h) {
     return x * x * u[0] + y * y * u[1] + z * z * u[2] + 2 * (x * y * u[3] + x * z * u[4] + y * z * u[5]);
 }
 
-// Each atom placed: its grid coordinates (fractional times the grid size, wrapped into the cell) for
-// the gather and, for a mask atom, its wrapped fractional coordinates and mask radius for the mask.
+// Each atom placed: its grid coordinates (fractional times the grid size, not wrapped, as the CPU's density
+// walk centres an atom) for the gather and, for a mask atom, its wrapped fractional coordinates and mask
+// radius for the mask.
 __global__ void place_kernel(const double *rel, int n, Placement p, int nu, int nv, int nw, const int *mask_slot,
                              const float *mask_radius, float4 *pos, double *mask_atoms) {
     const int i = blockIdx.x * blockDim.x + threadIdx.x;
@@ -96,20 +97,16 @@ __global__ void place_kernel(const double *rel, int n, Placement p, int nu, int
     const double c[3] = {p.r[0] * x + p.r[1] * y + p.r[2] * z + p.t[0],
                          p.r[3] * x + p.r[4] * y + p.r[5] * z + p.t[1],
                          p.r[6] * x + p.r[7] * y + p.r[8] * z + p.t[2]};
+    const int n3[3] = {nu, nv, nw};
     double f[3];
+    float g[3];
     for (int k = 0; k < 3; k++) {
         f[k] = p.frac[3 * k] * c[0] + p.frac[3 * k + 1] * c[1] + p.frac[3 * k + 2] * c[2];
+        g[k] = static_cast<float>(f[k] * n3[k]);
         f[k] -= floor(f[k]);
         if (f[k] >= 1.0)   // a tiny negative coordinate wraps to 1 - epsilon, which rounds to 1
             f[k] = 0.0;
     }
-    const int n3[3] = {nu, nv, nw};
-    float g[3];
-    for (int k = 0; k < 3; k++) {
-        g[k] = static_cast<float>(f[k] * n3[k]);
-        if (g[k] >= n3[k])
-            g[k] -= n3[k];
-    }
     pos[i] = make_float4(g[0], g[1], g[2], 0.0f);
     const int slot = mask_slot[i];
     if (slot >= 0) {

tests/CMakeLists.txt

@@ -154,6 +154,7 @@ ADD_EXECUTABLE(jfjoch_test
         ModelValidationTest.cpp
         MolecularReplacementTest.cpp
         ModelMaskGPUTest.cpp
+        ModelDensityGPUTest.cpp
         SearchSpaceGroupTest.cpp
         SearchSpaceGroupTwinTest.cpp
         SyntheticMergedReflections.h

tests/ModelDensityGPUTest.cpp

@@ -0,0 +1,122 @@
+// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
+// SPDX-License-Identifier: GPL-3.0-only
+
+#include <catch2/catch_all.hpp>
+#include "../common/CUDAWrapper.h"
+
+#ifdef JFJOCH_USE_CUDA
+
+#include <array>
+#include <cstring>
+#include <random>
+#include <vector>
+
+#include <gemmi/dencalc.hpp>
+#include <gemmi/it92.hpp>
+
+#include "../image_analysis/structure_refinement/ModelDensityGPU.h"
+#include "../image_analysis/structure_refinement/ModelGrid.h"
+#include "../image_analysis/structure_refinement/ModelStructureFactorsGPU.h"
+
+// The CUDA gather is the CPU's density walk, bit for bit (ModelDensityGPU.h): PutModelDensityOnGrid() with no
+// orbits to symmetrize - one copy of the content, which is what the gather puts on the grid - against
+// ModelDensityGPU fed as ModelStructureFactorsGPU feeds it. The grids are compared as bit patterns, as in
+// ModelDensityMetalTest.cpp.
+
+namespace {
+    // Atoms scattered over the whole cell and a half cell beyond it on every side - so that boxes wrap at every
+    // face and centres lie outside the cell - of a few elements, a quarter of them anisotropic, with occupancies
+    // below one.
+    gemmi::Structure random_crystal(const gemmi::UnitCell &cell, const char *group) {
+        std::mt19937 rng(11);
+        std::uniform_real_distribution<double> frac(-0.5, 1.5), b(8.0, 60.0), occ(0.3, 1.0), u(0.05, 0.4);
+        const gemmi::El elements[] = {gemmi::El::C, gemmi::El::N, gemmi::El::O, gemmi::El::S};
+        gemmi::Structure st;
+        st.cell = cell;
+        st.spacegroup_hm = group;
+        st.setup_cell_images();
+        gemmi::Model model(1);
+        gemmi::Chain chain("A");
+        gemmi::Residue res;
+        res.name = "UNK";
+        res.seqid = gemmi::SeqId(1, ' ');
+        for (int i = 0; i < 3000; i++) {
+            gemmi::Atom a;
+            a.name = "X";
+            a.element = elements[i % 4];
+            a.pos = cell.orthogonalize(gemmi::Fractional(frac(rng), frac(rng), frac(rng)));
+            a.b_iso = static_cast<float>(b(rng));
+            a.occ = static_cast<float>(occ(rng));
+            if (i % 4 == 3) {
+                const float d1 = static_cast<float>(u(rng)), d2 = static_cast<float>(u(rng)),
+                            d3 = static_cast<float>(u(rng));
+                a.aniso = {d1, d2, d3, 0.3f * std::min(d1, d2), -0.2f * std::min(d1, d3), 0.1f * std::min(d2, d3)};
+            }
+            res.atoms.push_back(a);
+        }
+        chain.residues.push_back(res);
+        model.chains.push_back(chain);
+        st.models.push_back(model);
+        return st;
+    }
+
+    void check_same_grid(const gemmi::UnitCell &cell, const char *group, double d_min) {
+        const gemmi::Structure st = random_crystal(cell, group);
+        gemmi::DensityCalculator<gemmi::IT92<float>, float> dc;
+        dc.d_min = d_min;
+        dc.rate = 1.5;
+        dc.set_grid_cell_and_spacegroup(st);
+        dc.set_refmac_compatible_blur(st.models[0]);
+        PutModelDensityOnGrid(dc, st.models[0], {}, 4, false);
+        const gemmi::Grid<float> &cpu = dc.grid;
+
+        ModelDensityGrid grid;
+        grid.nu = cpu.nu;
+        grid.nv = cpu.nv;
+        grid.nw = cpu.nw;
+        for (int i = 0; i < 3; i++)
+            for (int j = 0; j < 3; j++) {
+                grid.orth[3 * i + j] = cell.orth.mat[i][j];
+                grid.frac[3 * i + j] = cell.frac.mat[i][j];
+            }
+        std::vector<ModelDensityAtom> atoms;
+        std::vector<int> mask_atom;
+        std::vector<float> mask_radius;
+        ModelDensityAtoms(st.models[0], dc, atoms, mask_atom, mask_radius);
+        REQUIRE(ModelDensityGPU::Supports(grid, atoms));
+        // ModelStructureFactorsGPU::Compute()'s centres: the fractional position times the grid size.
+        std::vector<std::array<float, 4>> pos;
+        for (const gemmi::Atom &a : st.models[0].chains[0].residues[0].atoms) {
+            const gemmi::Fractional f = cell.fractionalize(a.pos);
+            pos.push_back({static_cast<float>(f.x * grid.nu), static_cast<float>(f.y * grid.nv),
+                           static_cast<float>(f.z * grid.nw), 0.0f});
+        }
+
+        CudaStream stream;
+        ModelDensityGPU density(stream, atoms.size(), ModelDensityGPU::PairBound(grid, atoms),
+                                ModelDensityGPU::Bricks(grid.nu, grid.nv, grid.nw));
+        density.SetAtoms(grid, atoms);
+        const size_t n = cpu.data.size();
+        CudaDevicePtr<float4> pos_d(pos.size());
+        CudaDevicePtr<float> grid_d(n);
+        REQUIRE(cudaMemcpyAsync(pos_d, pos.data(), pos.size() * sizeof(float4), cudaMemcpyHostToDevice, stream) == cudaSuccess);
+        density.Compute(pos_d, grid_d);
+        std::vector<float> gpu(n);
+        REQUIRE(cudaMemcpyAsync(gpu.data(), grid_d, n * sizeof(float), cudaMemcpyDeviceToHost, stream) == cudaSuccess);
+        REQUIRE(cudaStreamSynchronize(stream) == cudaSuccess);
+        size_t different = 0;
+        for (size_t i = 0; i < n; i++)
+            different += std::memcmp(&cpu.data[i], &gpu[i], sizeof(float)) != 0;
+        CHECK(different == 0);
+    }
+}
+
+TEST_CASE("ModelDensityGPU_BitIdenticalToCPU", "[ModelValidation][gpu]") {
+    if (get_gpu_count() == 0)
+        SKIP("No GPU");
+    // Triclinic, so that every term of the orthogonalization is there; then an orthorhombic group's grid.
+    check_same_grid(gemmi::UnitCell(61, 68, 77, 81, 96, 107), "P 1", 2.0);
+    check_same_grid(gemmi::UnitCell(70, 80, 90, 90, 90, 90), "P 21 21 21", 2.0);
+}
+
+#endif