src/gromacs/ewald/pme_spread.clh

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  *  \brief Implements PME OpenCL spline parameter computation and charge spread kernels.
  38  * When including this and other PME OpenCL kernel files, plenty of common
  39  * constants/macros are expected to be defined (such as "order" which is PME interpolation order).
  40  * For details, please see how pme_program.cl is compiled in pme_gpu_program_impl_ocl.cpp.
  41  *
  42  * This file's kernels specifically expect the following definitions:
  43  *
  44  * - atomsPerBlock which expresses how many atoms are processed by a single work group
  45  * - order which is a PME interpolation order
  46  * - computeSplines and spreadCharges must evaluate to either true or false to specify which
  47  * kernel falvor is being compiled
  48  * - wrapX and wrapY must evaluate to either true or false to specify whether the grid overlap
  49  * in dimension X/Y is to be used
  50  *
  51  *  \author Aleksei Iupinov <a.yupinov@gmail.com>
  52  */
  53
  54 #include "gromacs/gpu_utils/vectype_ops.clh"
  55
  56 #include "pme_gpu_calculate_splines.clh"
  57 #include "pme_gpu_types.h"
  58
  59 /*
  60  * This define affects the spline calculation behaviour in the kernel.
  61  * 0: a single GPU thread handles a single dimension of a single particle (calculating and storing
  62  * (order) spline values and derivatives). 1: (order) threads do redundant work on this same task,
  63  * each one stores only a single theta and single dtheta into global arrays. The only efficiency
  64  * difference is less global store operations, countered by more redundant spline computation.
  65  *
  66  * TODO: estimate if this should be a boolean parameter (and add it to the unit test if so).
  67  */
  68 #define PME_GPU_PARALLEL_SPLINE 0
  69
  70 #ifndef COMPILE_SPREAD_HELPERS_ONCE
  71 #    define COMPILE_SPREAD_HELPERS_ONCE
  72
  73 /*! \brief
  74  * General purpose function for loading atom-related data from global to shared memory.
  75  *
  76  * \param[out] sm_destination    Local memory array for output.
  77  * \param[in]  gm_source         Global memory array for input.
  78  * \param[in]  dataCountPerAtom  Number of data elements per single
  79  * atom (e.g. DIM for an rvec coordinates array).
  80  *
  81  */
  82 inline void pme_gpu_stage_atom_data(__local float* __restrict__ sm_destination,
  83                                     __global const float* __restrict__ gm_source,
  84                                     const int dataCountPerAtom)
  85 {
  86     assert((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0) + get_local_id(0)
  87            < MAX_INT);
  88     const int threadLocalIndex =
  89             (int)((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0)
  90                   + get_local_id(0));
  91     const int localIndex      = threadLocalIndex;
  92     const int globalIndexBase = (int)get_group_id(XX) * atomsPerBlock * dataCountPerAtom;
  93     const int globalIndex     = globalIndexBase + localIndex;
  94     if (localIndex < atomsPerBlock * dataCountPerAtom)
  95     {
  96         assert(isfinite(float(gm_source[globalIndex])));
  97         sm_destination[localIndex] = gm_source[globalIndex];
  98     }
  99 }
 100
 101 /*! \brief
 102  * PME GPU spline parameter and gridline indices calculation.
 103  * This corresponds to the CPU functions calc_interpolation_idx() and make_bsplines().
 104  * First stage of the whole kernel.
 105  *
 106  * \param[in]  kernelParams             Input PME GPU data in constant memory.
 107  * \param[in]  atomIndexOffset          Starting atom index for the execution block w.r.t. global memory.
 108  * \param[in]  sm_coordinates           Atom coordinates in the shared memory (rvec)
 109  * \param[in]  sm_coefficients          Atom charges/coefficients in the shared memory.
 110  * \param[out] sm_theta                 Atom spline values in the shared memory.
 111  * \param[out] sm_gridlineIndices       Atom gridline indices in the shared memory.
 112  * \param[out] sm_fractCoords           Atom fractional coordinates (rvec)
 113  * \param[out] gm_theta                 Atom spline parameter values
 114  * \param[out] gm_dtheta                Atom spline parameter derivatives
 115  * \param[out] gm_gridlineIndices       Same as \p sm_gridlineIndices but in global memory.
 116  * \param[in]  gm_fractShiftsTable      Atom fractional coordinates correction table
 117  * \param[in]  gm_gridlineIndicesTable  Table with atom gridline indices
 118  */
 119 gmx_opencl_inline void calculate_splines(const struct PmeOpenCLKernelParams kernelParams,
 120                                          const int                          atomIndexOffset,
 121                                          __local const float* __restrict__ sm_coordinates,
 122                                          __local const float* __restrict__ sm_coefficients,
 123                                          __local float* __restrict__ sm_theta,
 124                                          __local int* __restrict__ sm_gridlineIndices,
 125                                          __local float* __restrict__ sm_fractCoords,
 126                                          __global float* __restrict__ gm_theta,
 127                                          __global float* __restrict__ gm_dtheta,
 128                                          __global int* __restrict__ gm_gridlineIndices,
 129                                          __global const float* __restrict__ gm_fractShiftsTable,
 130                                          __global const int* __restrict__ gm_gridlineIndicesTable)
 131 {
 132     /* Thread index w.r.t. block */
 133     assert((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0) + get_local_id(0)
 134            < MAX_INT);
 135     const int threadLocalIndex =
 136             (int)((get_local_id(2) * get_local_size(1) + get_local_id(1)) * get_local_size(0)
 137                   + get_local_id(0));
 138     /* Warp index w.r.t. block - could probably be obtained easier? */
 139     const int warpIndex = threadLocalIndex / warp_size;
 140     /* Thread index w.r.t. warp */
 141     const int threadWarpIndex = threadLocalIndex % warp_size;
 142     /* Atom index w.r.t. warp - alternating 0 1 0 1 .. */
 143     const int atomWarpIndex = threadWarpIndex % atomsPerWarp;
 144     /* Atom index w.r.t. block/shared memory */
 145     const int atomIndexLocal = warpIndex * atomsPerWarp + atomWarpIndex;
 146
 147     /* Spline contribution index in one dimension */
 148     const int orderIndex = threadWarpIndex / (atomsPerWarp * DIM);
 149     /* Dimension index */
 150     const int dimIndex = (threadWarpIndex - orderIndex * (atomsPerWarp * DIM)) / atomsPerWarp;
 151
 152     /* Multi-purpose index of rvec/ivec atom data */
 153     const int sharedMemoryIndex = atomIndexLocal * DIM + dimIndex;
 154
 155     /* Spline parameters need a working buffer.
 156      * With PME_GPU_PARALLEL_SPLINE == 0 it is just a local array of (order) values for some of the threads, which is fine;
 157      * With PME_GPU_PARALLEL_SPLINE == 1 (order) times more threads are involved, so the shared memory is used to avoid overhead.
 158      * The buffer's size, striding and indexing are adjusted accordingly.
 159      * The buffer is accessed with SPLINE_DATA_PTR and SPLINE_DATA macros.
 160      */
 161 #    if PME_GPU_PARALLEL_SPLINE
 162 #        define splineDataStride (atomsPerBlock * DIM)
 163     const int      splineDataIndex = sharedMemoryIndex;
 164     __local float  sm_splineData[splineDataStride * order];
 165     __local float* splineDataPtr = sm_splineData;
 166 #    else
 167 #        define splineDataStride 1
 168     const int splineDataIndex = 0;
 169     float     splineData[splineDataStride * order];
 170     float*    splineDataPtr = splineData;
 171 #    endif
 172
 173 #    define SPLINE_DATA_PTR(i) (splineDataPtr + ((i)*splineDataStride + splineDataIndex))
 174 #    define SPLINE_DATA(i) (*SPLINE_DATA_PTR(i))
 175
 176     const int localCheck = (dimIndex < DIM) && (orderIndex < (PME_GPU_PARALLEL_SPLINE ? order : 1));
 177
 178     if (localCheck)
 179     {
 180         /* Indices interpolation */
 181         if (orderIndex == 0)
 182         {
 183             int          tableIndex = 0;
 184             float        n          = 0.0F;
 185             float        t          = 0.0F;
 186             const float3 x          = vload3(atomIndexLocal, sm_coordinates);
 187
 188             /* Accessing fields in fshOffset/nXYZ/recipbox/... with dimIndex offset
 189              * puts them into local memory(!) instead of accessing the constant memory directly.
 190              * That's the reason for the switch, to unroll explicitly.
 191              * The commented parts correspond to the 0 components of the recipbox.
 192              */
 193             switch (dimIndex)
 194             {
 195                 case XX:
 196                     tableIndex = kernelParams.grid.tablesOffsets[XX];
 197                     n          = kernelParams.grid.realGridSizeFP[XX];
 198                     t          = x.x * kernelParams.current.recipBox[dimIndex][XX]
 199                         + x.y * kernelParams.current.recipBox[dimIndex][YY]
 200                         + x.z * kernelParams.current.recipBox[dimIndex][ZZ];
 201                     break;
 202
 203                 case YY:
 204                     tableIndex = kernelParams.grid.tablesOffsets[YY];
 205                     n          = kernelParams.grid.realGridSizeFP[YY];
 206                     t          = /*x.x * kernelParams.current.recipbox[dimIndex][XX] + */ x.y
 207                                 * kernelParams.current.recipBox[dimIndex][YY]
 208                         + x.z * kernelParams.current.recipBox[dimIndex][ZZ];
 209                     break;
 210
 211                 case ZZ:
 212                     tableIndex = kernelParams.grid.tablesOffsets[ZZ];
 213                     n          = kernelParams.grid.realGridSizeFP[ZZ];
 214                     t          = /*x.x * kernelParams.current.recipbox[dimIndex][XX] + x.y * kernelParams.current.recipbox[dimIndex][YY] + */ x
 215                                 .z
 216                         * kernelParams.current.recipBox[dimIndex][ZZ];
 217                     break;
 218                 default:
 219                     assert(false);
 220                     return;
 221                     break;
 222             }
 223             const float shift = c_pmeMaxUnitcellShift;
 224
 225             /* Fractional coordinates along box vectors, adding a positive shift to ensure t is positive for triclinic boxes */
 226             t                                 = (t + shift) * n;
 227             const int tInt                    = (int)t;
 228             sm_fractCoords[sharedMemoryIndex] = t - (float)tInt;
 229             tableIndex += tInt;
 230             assert(tInt >= 0);
 231             assert(tInt < c_pmeNeighborUnitcellCount * n);
 232             sm_fractCoords[sharedMemoryIndex] += gm_fractShiftsTable[tableIndex];
 233             sm_gridlineIndices[sharedMemoryIndex] = gm_gridlineIndicesTable[tableIndex];
 234             gm_gridlineIndices[atomIndexOffset * DIM + sharedMemoryIndex] =
 235                     sm_gridlineIndices[sharedMemoryIndex];
 236         }
 237
 238         /* B-spline calculation */
 239
 240         const int chargeCheck = pme_gpu_check_atom_charge(sm_coefficients[atomIndexLocal]);
 241         if (chargeCheck)
 242         {
 243             int o = orderIndex; // This is an index that is set once for PME_GPU_PARALLEL_SPLINE == 1
 244
 245             const float dr = sm_fractCoords[sharedMemoryIndex];
 246             assert(isfinite(dr));
 247
 248             /* dr is relative offset from lower cell limit */
 249             *SPLINE_DATA_PTR(order - 1) = 0.0F;
 250             *SPLINE_DATA_PTR(1)         = dr;
 251             *SPLINE_DATA_PTR(0)         = 1.0F - dr;
 252
 253 #    pragma unroll order
 254             for (int k = 3; k < order; k++)
 255             {
 256                 const float div         = 1.0F / ((float)k - 1.0F);
 257                 *SPLINE_DATA_PTR(k - 1) = div * dr * SPLINE_DATA(k - 2);
 258 #    pragma unroll
 259                 for (int l = 1; l < (k - 1); l++)
 260                 {
 261                     *SPLINE_DATA_PTR(k - l - 1) =
 262                             div
 263                             * ((dr + (float)l) * SPLINE_DATA(k - l - 2)
 264                                + ((float)k - (float)l - dr) * SPLINE_DATA(k - l - 1));
 265                 }
 266                 *SPLINE_DATA_PTR(0) = div * (1.0F - dr) * SPLINE_DATA(0);
 267             }
 268
 269             const int thetaIndexBase        = getSplineParamIndexBase(warpIndex, atomWarpIndex);
 270             const int thetaGlobalOffsetBase = atomIndexOffset * DIM * order;
 271
 272             /* Differentiation and storing the spline derivatives (dtheta) */
 273 #    if !PME_GPU_PARALLEL_SPLINE
 274             // With PME_GPU_PARALLEL_SPLINE == 1, o is already set to orderIndex;
 275             // With PME_GPU_PARALLEL_SPLINE == 0, we loop o over range(order).
 276 #        pragma unroll order
 277             for (o = 0; o < order; o++)
 278 #    endif
 279             {
 280                 const int thetaIndex       = getSplineParamIndex(thetaIndexBase, dimIndex, o);
 281                 const int thetaGlobalIndex = thetaGlobalOffsetBase + thetaIndex;
 282
 283                 const float dtheta = ((o > 0) ? SPLINE_DATA(o - 1) : 0.0F) - SPLINE_DATA(o);
 284                 assert(isfinite(dtheta));
 285                 gm_dtheta[thetaGlobalIndex] = dtheta;
 286             }
 287
 288             const float div             = 1.0F / (order - 1.0F);
 289             *SPLINE_DATA_PTR(order - 1) = div * dr * SPLINE_DATA(order - 2);
 290 #    pragma unroll
 291             for (int k = 1; k < (order - 1); k++)
 292             {
 293                 *SPLINE_DATA_PTR(order - k - 1) = div
 294                                                   * ((dr + (float)k) * SPLINE_DATA(order - k - 2)
 295                                                      + (order - k - dr) * SPLINE_DATA(order - k - 1));
 296             }
 297             *SPLINE_DATA_PTR(0) = div * (1.0F - dr) * SPLINE_DATA(0);
 298
 299             /* Storing the spline values (theta) */
 300 #    if !PME_GPU_PARALLEL_SPLINE
 301             // See comment for the loop above
 302 #        pragma unroll order
 303             for (o = 0; o < order; o++)
 304 #    endif
 305             {
 306                 const int thetaIndex       = getSplineParamIndex(thetaIndexBase, dimIndex, o);
 307                 const int thetaGlobalIndex = thetaGlobalOffsetBase + thetaIndex;
 308
 309                 sm_theta[thetaIndex] = SPLINE_DATA(o);
 310                 assert(isfinite(sm_theta[thetaIndex]));
 311                 gm_theta[thetaGlobalIndex] = SPLINE_DATA(o);
 312             }
 313         }
 314     }
 315 }
 316
 317 /*! \brief
 318  * Charge spreading onto the grid.
 319  * This corresponds to the CPU function spread_coefficients_bsplines_thread().
 320  * Second stage of the whole kernel.
 321  *
 322  * \param[in]  kernelParams         Input PME GPU data in constant memory.
 323  * \param[in]  sm_coefficients      Atom coefficents/charges in the shared memory.
 324  * \param[in]  sm_gridlineIndices   Atom gridline indices in the shared memory.
 325  * \param[in]  sm_theta             Atom spline values in the shared memory.
 326  * \param[out] gm_grid              Global 3D grid for spreading.
 327  */
 328 gmx_opencl_inline void spread_charges(const struct PmeOpenCLKernelParams kernelParams,
 329                                       __local const float* __restrict__ sm_coefficients,
 330                                       __local const int* __restrict__ sm_gridlineIndices,
 331                                       __local const float* __restrict__ sm_theta,
 332                                       __global float* __restrict__ gm_grid)
 333 {
 334     const int nx  = kernelParams.grid.realGridSize[XX];
 335     const int ny  = kernelParams.grid.realGridSize[YY];
 336     const int nz  = kernelParams.grid.realGridSize[ZZ];
 337     const int pny = kernelParams.grid.realGridSizePadded[YY];
 338     const int pnz = kernelParams.grid.realGridSizePadded[ZZ];
 339
 340     const int offx = 0;
 341     const int offy = 0;
 342     const int offz = 0;
 343
 344     const int atomIndexLocal = get_local_id(ZZ);
 345
 346     const int chargeCheck = pme_gpu_check_atom_charge(sm_coefficients[atomIndexLocal]);
 347     if (chargeCheck)
 348     {
 349         // Spline Y/Z coordinates
 350         const int ithy   = get_local_id(YY);
 351         const int ithz   = get_local_id(XX);
 352         const int ixBase = sm_gridlineIndices[atomIndexLocal * DIM + XX] - offx;
 353         int       iy     = sm_gridlineIndices[atomIndexLocal * DIM + YY] - offy + ithy;
 354         if (wrapY & (iy >= ny))
 355         {
 356             iy -= ny;
 357         }
 358         int iz = sm_gridlineIndices[atomIndexLocal * DIM + ZZ] - offz + ithz;
 359         if (iz >= nz)
 360         {
 361             iz -= nz;
 362         }
 363
 364         /* Atom index w.r.t. warp - alternating 0 1 0 1 .. */
 365         const int atomWarpIndex = atomIndexLocal % atomsPerWarp;
 366         /* Warp index w.r.t. block - could probably be obtained easier? */
 367         const int warpIndex = atomIndexLocal / atomsPerWarp;
 368
 369         const int   splineIndexBase = getSplineParamIndexBase(warpIndex, atomWarpIndex);
 370         const int   splineIndexZ    = getSplineParamIndex(splineIndexBase, ZZ, ithz);
 371         const float thetaZ          = sm_theta[splineIndexZ];
 372         const int   splineIndexY    = getSplineParamIndex(splineIndexBase, YY, ithy);
 373         const float thetaY          = sm_theta[splineIndexY];
 374         const float constVal        = thetaZ * thetaY * sm_coefficients[atomIndexLocal];
 375         assert(isfinite(constVal));
 376         const int constOffset = iy * pnz + iz;
 377
 378 #    pragma unroll order
 379         for (int ithx = 0; (ithx < order); ithx++)
 380         {
 381             int ix = ixBase + ithx;
 382             if (wrapX & (ix >= nx))
 383             {
 384                 ix -= nx;
 385             }
 386             const int   gridIndexGlobal = ix * pny * pnz + constOffset;
 387             const int   splineIndexX    = getSplineParamIndex(splineIndexBase, XX, ithx);
 388             const float thetaX          = sm_theta[splineIndexX];
 389             assert(isfinite(thetaX));
 390             assert(isfinite(gm_grid[gridIndexGlobal]));
 391             atomicAdd_g_f(gm_grid + gridIndexGlobal, thetaX * constVal);
 392         }
 393     }
 394 }
 395
 396 #endif // COMPILE_SPREAD_HELPERS_ONCE
 397
 398 /*! \brief
 399  * A spline computation and/or charge spreading kernel function.
 400  * Please see the file description for additional defines which this kernel expects.
 401  *
 402  * \param[in]     kernelParams             Input PME GPU data in constant memory.
 403  * \param[in,out] gm_theta                 Atom spline parameter values
 404  * \param[out]    gm_dtheta                Atom spline parameter derivatives
 405  * \param[in,out] gm_gridlineIndices       Atom gridline indices (ivec)
 406  * \param[out]    gm_gridA                 Global 3D grid for charge spreading.
 407  * Grid for the unperturbed state, FEP state A or the single grid used for
 408  * interpolated coefficients on one grid in FEP A/B.
 409  * \param[out]    gm_gridB                 Global 3D grid for charge spreading.
 410  * FEP state B when using dual grids (when calculating energy and virials).
 411  * \param[in]     gm_fractShiftsTable      Atom fractional coordinates correction table
 412  * \param[in]     gm_gridlineIndicesTable  Table with atom gridline indices
 413  * \param[in]     gm_coefficientsA         Atom charges/coefficients of the unperturbed state
 414  * or FEP state A.
 415  * \param[in]     gm_coefficientsB         Atom charges/coefficients in FEP state B. Only
 416  * used when spreading interpolated coefficients on one grid.
 417  * \param[in]     gm_coordinates           Atom coordinates (rvec)
 418  */
 419 __attribute__((reqd_work_group_size(order, order, atomsPerBlock))) __kernel void CUSTOMIZED_KERNEL_NAME(
 420         pme_spline_and_spread_kernel)(const struct PmeOpenCLKernelParams kernelParams,
 421                                       __global float* __restrict__ gm_theta,
 422                                       __global float* __restrict__ gm_dtheta,
 423                                       __global int* __restrict__ gm_gridlineIndices,
 424                                       __global float* __restrict__ gm_gridA,
 425                                       __global float* __restrict__ gm_gridB,
 426                                       __global const float* __restrict__ gm_fractShiftsTable,
 427                                       __global const int* __restrict__ gm_gridlineIndicesTable,
 428                                       __global const float* __restrict__ gm_coefficientsA,
 429                                       __global const float* __restrict__ gm_coefficientsB,
 430                                       __global const float* __restrict__ gm_coordinates)
 431 {
 432     // Gridline indices, ivec
 433     __local int sm_gridlineIndices[atomsPerBlock * DIM];
 434     // Charges
 435     __local float sm_coefficients[atomsPerBlock];
 436     // Spline values
 437     __local float sm_theta[atomsPerBlock * DIM * order];
 438     // Fractional coordinates - only for spline computation
 439     __local float sm_fractCoords[atomsPerBlock * DIM];
 440     // Staging coordinates - only for spline computation
 441     __local float sm_coordinates[DIM * atomsPerBlock];
 442
 443     const int atomIndexOffset = (int)get_group_id(XX) * atomsPerBlock;
 444     int       gridIndex       = 0;
 445
 446     /* Staging coefficients/charges for both spline and spread */
 447     pme_gpu_stage_atom_data(sm_coefficients, gm_coefficientsA, 1);
 448
 449     if (computeSplines)
 450     {
 451         /* Staging coordinates */
 452         pme_gpu_stage_atom_data(sm_coordinates, gm_coordinates, DIM);
 453
 454         barrier(CLK_LOCAL_MEM_FENCE);
 455         calculate_splines(kernelParams, atomIndexOffset, sm_coordinates, sm_coefficients, sm_theta,
 456                           sm_gridlineIndices, sm_fractCoords, gm_theta, gm_dtheta,
 457                           gm_gridlineIndices, gm_fractShiftsTable, gm_gridlineIndicesTable);
 458 #if !defined(_AMD_SOURCE_) && !defined(_NVIDIA_SOURCE_)
 459         /* This is only here for execution of e.g. 32-sized warps on 16-wide hardware; this was
 460          * __syncwarp() in CUDA. #2519
 461          */
 462         barrier(CLK_LOCAL_MEM_FENCE);
 463 #endif
 464     }
 465     else
 466     {
 467         /* Staging the data for spread
 468          * (the data is assumed to be in GPU global memory with proper layout already,
 469          * as in after running the spline kernel)
 470          */
 471         /* Spline data - only thetas (dthetas will only be needed in gather) */
 472         pme_gpu_stage_atom_data(sm_theta, gm_theta, DIM * order);
 473         /* Gridline indices - they're actually int and not float, but C99 is angry about overloads */
 474         pme_gpu_stage_atom_data((__local float*)sm_gridlineIndices,
 475                                 (__global const float*)gm_gridlineIndices, DIM);
 476
 477         barrier(CLK_LOCAL_MEM_FENCE);
 478     }
 479     /* Spreading */
 480     if (spreadCharges)
 481     {
 482         spread_charges(kernelParams, sm_coefficients, sm_gridlineIndices, sm_theta, gm_gridA);
 483     }
 484     if (numGrids == 2)
 485     {
 486         barrier(CLK_LOCAL_MEM_FENCE);
 487         gridIndex = 1;
 488         pme_gpu_stage_atom_data(sm_coefficients, gm_coefficientsB, 1);
 489         barrier(CLK_LOCAL_MEM_FENCE);
 490         if (spreadCharges)
 491         {
 492             spread_charges(kernelParams, sm_coefficients, sm_gridlineIndices, sm_theta, gm_gridB);
 493         }
 494     }
 495 }