diff --git a/image_analysis/structure_refinement/ModelDensityGPU.cu b/image_analysis/structure_refinement/ModelDensityGPU.cu index d85f9e642..30bd7dd38 100644 --- a/image_analysis/structure_refinement/ModelDensityGPU.cu +++ b/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(b.lo[0] + mu); + const float fy = b.nctr[1] - static_cast(b.lo[1] + mv); + const float fz = b.nctr[2] - static_cast(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 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(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(grid.orth[0] * inv_n[0]); + geom.m01 = static_cast(grid.orth[1] * inv_n[1]); + geom.m02 = static_cast(grid.orth[2] * inv_n[2]); + geom.m11 = static_cast(grid.orth[4] * inv_n[1]); + geom.m12 = static_cast(grid.orth[5] * inv_n[2]); + geom.m22 = static_cast(grid.orth[8] * inv_n[2]); const std::vector 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(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<<>>(key_sorted, npairs, brick_start, brick_end); cuda_err(cudaGetLastError()); } - gather_kernel<<>>(atoms, d_pos, value_sorted, brick_start, brick_end, + gather_kernel<<>>(atoms, d_pos, box, value_sorted, brick_start, brick_end, geom, d_grid); cuda_err(cudaGetLastError()); } diff --git a/image_analysis/structure_refinement/ModelDensityGPU.h b/image_analysis/structure_refinement/ModelDensityGPU.h index d58734d9c..d46941ab7 100644 --- a/image_analysis/structure_refinement/ModelDensityGPU.h +++ b/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 #include @@ -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 &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: diff --git a/image_analysis/structure_refinement/ModelStructureFactorsGPU.cpp b/image_analysis/structure_refinement/ModelStructureFactorsGPU.cpp index a86b8eb9a..7f2f9373a 100644 --- a/image_analysis/structure_refinement/ModelStructureFactorsGPU.cpp +++ b/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(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(radius); + a.max_r2 = static_cast(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(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(radius); + a.max_r2 = static_cast(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 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> pos; std::vector> frac; @@ -202,12 +207,10 @@ void ModelStructureFactorsGPU::Compute(const gemmi::Model &model, gemmi::AsuData std::array f{f0.x, f0.y, f0.z}; std::array g{}; for (int k = 0; k < 3; k++) { + g[k] = static_cast(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(f[k] * n3[k]); - if (g[k] >= n3[k]) - g[k] -= n3[k]; } frac.push_back(f); pos.push_back(g); diff --git a/image_analysis/structure_refinement/RigidBodyGPU.cu b/image_analysis/structure_refinement/RigidBodyGPU.cu index 2b73e5b3d..962221384 100644 --- a/image_analysis/structure_refinement/RigidBodyGPU.cu +++ b/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(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(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) { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 1895e035d..8a670baf0 100644 --- a/tests/CMakeLists.txt +++ b/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 diff --git a/tests/ModelDensityGPUTest.cpp b/tests/ModelDensityGPUTest.cpp new file mode 100644 index 000000000..fac8cec33 --- /dev/null +++ b/tests/ModelDensityGPUTest.cpp @@ -0,0 +1,122 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "../common/CUDAWrapper.h" + +#ifdef JFJOCH_USE_CUDA + +#include +#include +#include +#include + +#include +#include + +#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 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(b(rng)); + a.occ = static_cast(occ(rng)); + if (i % 4 == 3) { + const float d1 = static_cast(u(rng)), d2 = static_cast(u(rng)), + d3 = static_cast(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, 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 &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 atoms; + std::vector mask_atom; + std::vector 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> 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(f.x * grid.nu), static_cast(f.y * grid.nv), + static_cast(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 pos_d(pos.size()); + CudaDevicePtr grid_d(n); + REQUIRE(cudaMemcpyAsync(pos_d, pos.data(), pos.size() * sizeof(float4), cudaMemcpyHostToDevice, stream) == cudaSuccess); + density.Compute(pos_d, grid_d); + std::vector 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