src/gromacs/nbnxm/opencl/nbnxm_ocl_kernel_utils.clh

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2016,2017 by the GROMACS development team.
   5  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  *  \brief
  38  *  Utility constant and function declaration for the OpenCL non-bonded kernels.
  39  *  This header should be included once at the top level, just before the
  40  *  kernels are included (has to be preceded by nbnxn_ocl_types.h).
  41  *
  42  *  \author Szilárd Páll <pall.szilard@gmail.com>
  43  *  \ingroup module_nbnxm
  44  */
  45
  46 #define GMX_DOUBLE 0
  47
  48 #include "gromacs/gpu_utils/device_utils.clh"
  49 #include "gromacs/gpu_utils/vectype_ops.clh"
  50 #include "gromacs/pbcutil/ishift.h"
  51
  52 #include "nbnxm_ocl_consts.h"
  53
  54 #define CL_SIZE (c_nbnxnGpuClusterSize)
  55
  56 #define WARP_SIZE (CL_SIZE * CL_SIZE / 2) // Currently only c_nbnxnGpuClusterpairSplit=2 supported
  57
  58 #if defined _NVIDIA_SOURCE_ || defined _AMD_SOURCE_
  59 /* Currently we enable CJ prefetch for AMD/NVIDIA and disable it for other vendors
  60  * Note that this should precede the kernel_utils include.
  61  */
  62 #    define USE_CJ_PREFETCH 1
  63 #else
  64 #    define USE_CJ_PREFETCH 0
  65 #endif
  66
  67 #if defined cl_intel_subgroups || defined cl_khr_subgroups \
  68         || (defined __OPENCL_VERSION__ && __OPENCL_VERSION__ >= 210)
  69 #    define HAVE_SUBGROUP 1
  70 #else
  71 #    define HAVE_SUBGROUP 0
  72 #endif
  73
  74 #ifdef cl_intel_subgroups
  75 #    define HAVE_INTEL_SUBGROUP 1
  76 #else
  77 #    define HAVE_INTEL_SUBGROUP 0
  78 #endif
  79
  80 #if defined _INTEL_SOURCE_
  81 #    define SUBGROUP_SIZE 8
  82 #elif defined _AMD_SOURCE_
  83 #    define SUBGROUP_SIZE 64
  84 #else
  85 #    define SUBGROUP_SIZE 32
  86 #endif
  87
  88 #define REDUCE_SHUFFLE (HAVE_INTEL_SUBGROUP && CL_SIZE == 4 && SUBGROUP_SIZE == WARP_SIZE)
  89 #define USE_SUBGROUP_ANY (HAVE_SUBGROUP && SUBGROUP_SIZE == WARP_SIZE)
  90 #define USE_SUBGROUP_PRELOAD HAVE_INTEL_SUBGROUP
  91
  92 /* 1.0 / sqrt(M_PI) */
  93 #define M_FLOAT_1_SQRTPI 0.564189583547756F
  94
  95 //-------------------
  96
  97 #ifndef NBNXN_OPENCL_KERNEL_UTILS_CLH
  98 #    define NBNXN_OPENCL_KERNEL_UTILS_CLH
  99
 100 #    if CL_SIZE == 8
 101 #        define WARP_SIZE_LOG2 (5)
 102 #        define CL_SIZE_LOG2 (3)
 103 #    elif CL_SIZE == 4
 104 #        define WARP_SIZE_LOG2 (3)
 105 #        define CL_SIZE_LOG2 (2)
 106 #    else
 107 #        error unsupported CL_SIZE
 108 #    endif
 109
 110 #    define CL_SIZE_SQ (CL_SIZE * CL_SIZE)
 111 #    define FBUF_STRIDE (CL_SIZE_SQ)
 112
 113 #    define ONE_SIXTH_F 0.16666667F
 114 #    define ONE_TWELVETH_F 0.08333333F
 115 #    define HALF_F 0.5F
 116
 117
 118 #    ifdef __GNUC__
 119 /* GCC, clang, and some ICC pretending to be GCC */
 120 #        define gmx_unused __attribute__((unused))
 121 #    else
 122 #        define gmx_unused
 123 #    endif
 124
 125 // Data structures shared between OpenCL device code and OpenCL host code
 126 // TODO: review, improve
 127 // Replaced real by float for now, to avoid including any other header
 128 typedef struct
 129 {
 130     float c2;
 131     float c3;
 132     float cpot;
 133 } shift_consts_t;
 134
 135 /* Used with potential switching:
 136  * rsw        = max(r - r_switch, 0)
 137  * sw         = 1 + c3*rsw^3 + c4*rsw^4 + c5*rsw^5
 138  * dsw        = 3*c3*rsw^2 + 4*c4*rsw^3 + 5*c5*rsw^4
 139  * force      = force*dsw - potential*sw
 140  * potential *= sw
 141  */
 142 typedef struct
 143 {
 144     float c3;
 145     float c4;
 146     float c5;
 147 } switch_consts_t;
 148
 149 // Data structure shared between the OpenCL device code and OpenCL host code
 150 // Must not contain OpenCL objects (buffers)
 151 typedef struct cl_nbparam_params
 152 {
 153
 154     //! type of electrostatics, takes values from #eelCu
 155     int eeltype;
 156     //! type of VdW impl., takes values from #evdwCu
 157     int vdwtype;
 158
 159     //! charge multiplication factor
 160     float epsfac;
 161     //! Reaction-field/plain cutoff electrostatics const.
 162     float c_rf;
 163     //! Reaction-field electrostatics constant
 164     float two_k_rf;
 165     //! Ewald/PME parameter
 166     float ewald_beta;
 167     //! Ewald/PME correction term substracted from the direct-space potential
 168     float sh_ewald;
 169     //! LJ-Ewald/PME correction term added to the correction potential
 170     float sh_lj_ewald;
 171     //! LJ-Ewald/PME coefficient
 172     float ewaldcoeff_lj;
 173
 174     //! Coulomb cut-off squared
 175     float rcoulomb_sq;
 176
 177     //! VdW cut-off squared
 178     float rvdw_sq;
 179     //! VdW switched cut-off
 180     float rvdw_switch;
 181     //! Full, outer pair-list cut-off squared
 182     float rlistOuter_sq;
 183     //! Inner, dynamic pruned pair-list cut-off squared  XXX: this is only needed in the pruning kernels, but for now we also pass it to the nonbondeds
 184     float rlistInner_sq;
 185
 186     //! VdW shift dispersion constants
 187     shift_consts_t dispersion_shift;
 188     //! VdW shift repulsion constants
 189     shift_consts_t repulsion_shift;
 190     //! VdW switch constants
 191     switch_consts_t vdw_switch;
 192
 193     /* Ewald Coulomb force table data - accessed through texture memory */
 194     //! table scale/spacing
 195     float coulomb_tab_scale;
 196 } cl_nbparam_params_t;
 197
 198 typedef struct
 199 {
 200     //! i-super-cluster
 201     int sci;
 202     //! Shift vector index plus possible flags
 203     int shift;
 204     //! Start index into cj4
 205     int cj4_ind_start;
 206     //! End index into cj4
 207     int cj4_ind_end;
 208 } nbnxn_sci_t;
 209
 210 typedef struct
 211 {
 212     //! The i-cluster interactions mask for 1 warp
 213     unsigned int imask;
 214     //! Index into the exclusion array for 1 warp
 215     int excl_ind;
 216 } nbnxn_im_ei_t;
 217
 218 typedef struct
 219 {
 220     //! The 4 j-clusters
 221     int cj[4];
 222     //! The i-cluster mask data       for 2 warps
 223     nbnxn_im_ei_t imei[2];
 224 } nbnxn_cj4_t;
 225
 226
 227 typedef struct
 228 {
 229     unsigned int pair[CL_SIZE * CL_SIZE / 2]; /* Topology exclusion interaction bits for one warp,
 230                                                * each unsigned has bitS for 4*8 i clusters
 231                                                */
 232 } nbnxn_excl_t;
 233
 234 /*! i-cluster interaction mask for a super-cluster with all c_nbnxnGpuNumClusterPerSupercluster bits set */
 235 __constant unsigned supercl_interaction_mask = ((1U << c_nbnxnGpuNumClusterPerSupercluster) - 1U);
 236
 237 /*! Minimum single precision threshold for r^2 to avoid r^-12 overflow. */
 238 __constant float c_nbnxnMinDistanceSquared = NBNXM_MIN_DISTANCE_SQUARED_VALUE_FLOAT;
 239
 240 gmx_opencl_inline void preloadCj4Generic(__local int*        sm_cjPreload,
 241                                          const __global int* gm_cj,
 242                                          int                 tidxi,
 243                                          int                 tidxj,
 244                                          bool gmx_unused iMaskCond)
 245 {
 246     /* Pre-load cj into shared memory */
 247 #    if defined _AMD_SOURCE_ // TODO: fix by setting c_nbnxnGpuClusterpairSplit properly
 248     if (tidxj == 0 & tidxi < c_nbnxnGpuJgroupSize)
 249     {
 250         sm_cjPreload[tidxi] = gm_cj[tidxi];
 251     }
 252 #    else
 253     const int c_clSize                   = CL_SIZE;
 254     const int c_nbnxnGpuClusterpairSplit = 2;
 255     const int c_splitClSize              = c_clSize / c_nbnxnGpuClusterpairSplit;
 256     if ((tidxj == 0 | tidxj == c_splitClSize) & (tidxi < c_nbnxnGpuJgroupSize))
 257     {
 258         sm_cjPreload[tidxi + tidxj * c_nbnxnGpuJgroupSize / c_splitClSize] = gm_cj[tidxi];
 259     }
 260 #    endif
 261 }
 262
 263
 264 #    if USE_SUBGROUP_PRELOAD
 265 gmx_opencl_inline int preloadCj4Subgroup(const __global int* gm_cj)
 266 {
 267     // loads subgroup-size # of elements (8) instead of the 4 required
 268     // equivalent to *cjs = *gm_cj
 269     return intel_sub_group_block_read((const __global uint*)gm_cj);
 270 }
 271 #    endif // USE_SUBGROUP_PRELOAD
 272
 273 #    if USE_SUBGROUP_PRELOAD
 274 typedef size_t CjType;
 275 #    else
 276 typedef __local int* CjType;
 277 #    endif
 278
 279 /*! \brief Preload cj4
 280  *
 281  * - For AMD we load once for a wavefront of 64 threads (on 4 threads * NTHREAD_Z)
 282  * - For NVIDIA once per warp (on 2x4 threads * NTHREAD_Z)
 283  * - For Intel(/USE_SUBGROUP_PRELOAD) loads into private memory(/register) instead of local memory
 284  *
 285  * It is the caller's responsibility to make sure that data is consumed only when
 286  * it's ready. This function does not call a barrier.
 287  */
 288 gmx_opencl_inline void preloadCj4(CjType gmx_unused* cjs,
 289                                   const __global int gmx_unused* gm_cj,
 290                                   int gmx_unused tidxi,
 291                                   int gmx_unused tidxj,
 292                                   bool gmx_unused iMaskCond)
 293 {
 294 #    if USE_SUBGROUP_PRELOAD
 295     *cjs = preloadCj4Subgroup(gm_cj);
 296 #    elif USE_CJ_PREFETCH
 297     preloadCj4Generic(*cjs, gm_cj, tidxi, tidxj, iMaskCond);
 298 #    else
 299     // nothing to do
 300 #    endif
 301 }
 302
 303 gmx_opencl_inline int loadCjPreload(__local int* sm_cjPreload, int jm, int gmx_unused tidxi, int gmx_unused tidxj)
 304 {
 305 #    if defined _AMD_SOURCE_
 306     int         warpLoadOffset = 0; // TODO: fix by setting c_nbnxnGpuClusterpairSplit properly
 307 #    else
 308     const int c_clSize = CL_SIZE;
 309     const int c_nbnxnGpuClusterpairSplit = 2;
 310     const int c_splitClSize = c_clSize / c_nbnxnGpuClusterpairSplit;
 311     int warpLoadOffset = (tidxj & c_splitClSize) * c_nbnxnGpuJgroupSize / c_splitClSize;
 312 #    endif
 313     return sm_cjPreload[jm + warpLoadOffset];
 314 }
 315
 316 /* \brief Load a cj given a jm index.
 317  *
 318  * If cj4 preloading is enabled, it loads from the local memory, otherwise from global.
 319  */
 320 gmx_opencl_inline int
 321 loadCj(CjType cjs, const __global int gmx_unused* gm_cj, int jm, int gmx_unused tidxi, int gmx_unused tidxj)
 322 {
 323 #    if USE_SUBGROUP_PRELOAD
 324     return sub_group_broadcast(cjs, jm);
 325 #    elif USE_CJ_PREFETCH
 326     return loadCjPreload(cjs, jm, tidxi, tidxj);
 327 #    else
 328     return gm_cj[jm];
 329 #    endif
 330 }
 331
 332 /*! Convert LJ sigma,epsilon parameters to C6,C12. */
 333 gmx_opencl_inline float2 convert_sigma_epsilon_to_c6_c12(const float sigma, const float epsilon)
 334 {
 335     const float  sigma2 = sigma * sigma;
 336     const float  sigma6 = sigma2 * sigma2 * sigma2;
 337     const float  c6     = epsilon * sigma6;
 338     const float2 c6c12  = (float2)(c6,           /* c6 */
 339                                   c6 * sigma6); /* c12 */
 340     return c6c12;
 341 }
 342
 343
 344 /*! Apply force switch,  force + energy version. */
 345 gmx_opencl_inline void calculate_force_switch_F(const cl_nbparam_params_t* nbparam,
 346                                                 float                      c6,
 347                                                 float                      c12,
 348                                                 float                      inv_r,
 349                                                 float                      r2,
 350                                                 float*                     F_invr)
 351 {
 352     /* force switch constants */
 353     const float disp_shift_V2 = nbparam->dispersion_shift.c2;
 354     const float disp_shift_V3 = nbparam->dispersion_shift.c3;
 355     const float repu_shift_V2 = nbparam->repulsion_shift.c2;
 356     const float repu_shift_V3 = nbparam->repulsion_shift.c3;
 357
 358     const float r        = r2 * inv_r;
 359     float       r_switch = r - nbparam->rvdw_switch;
 360     r_switch             = r_switch >= 0.0F ? r_switch : 0.0F;
 361
 362     *F_invr += -c6 * (disp_shift_V2 + disp_shift_V3 * r_switch) * r_switch * r_switch * inv_r
 363                + c12 * (repu_shift_V2 + repu_shift_V3 * r_switch) * r_switch * r_switch * inv_r;
 364 }
 365
 366 /*! Apply force switch, force-only version. */
 367 gmx_opencl_inline void calculate_force_switch_F_E(const cl_nbparam_params_t* nbparam,
 368                                                   float                      c6,
 369                                                   float                      c12,
 370                                                   float                      inv_r,
 371                                                   float                      r2,
 372                                                   float*                     F_invr,
 373                                                   float*                     E_lj)
 374 {
 375     /* force switch constants */
 376     const float disp_shift_V2 = nbparam->dispersion_shift.c2;
 377     const float disp_shift_V3 = nbparam->dispersion_shift.c3;
 378     const float repu_shift_V2 = nbparam->repulsion_shift.c2;
 379     const float repu_shift_V3 = nbparam->repulsion_shift.c3;
 380
 381     const float disp_shift_F2 = nbparam->dispersion_shift.c2 / 3;
 382     const float disp_shift_F3 = nbparam->dispersion_shift.c3 / 4;
 383     const float repu_shift_F2 = nbparam->repulsion_shift.c2 / 3;
 384     const float repu_shift_F3 = nbparam->repulsion_shift.c3 / 4;
 385
 386     const float r        = r2 * inv_r;
 387     float       r_switch = r - nbparam->rvdw_switch;
 388     r_switch             = r_switch >= 0.0F ? r_switch : 0.0F;
 389
 390     *F_invr += -c6 * (disp_shift_V2 + disp_shift_V3 * r_switch) * r_switch * r_switch * inv_r
 391                + c12 * (repu_shift_V2 + repu_shift_V3 * r_switch) * r_switch * r_switch * inv_r;
 392     *E_lj += c6 * (disp_shift_F2 + disp_shift_F3 * r_switch) * r_switch * r_switch * r_switch
 393              - c12 * (repu_shift_F2 + repu_shift_F3 * r_switch) * r_switch * r_switch * r_switch;
 394 }
 395
 396 /*! Apply potential switch, force-only version. */
 397 gmx_opencl_inline void calculate_potential_switch_F(const cl_nbparam_params_t* nbparam,
 398                                                     float                      inv_r,
 399                                                     float                      r2,
 400                                                     float*                     F_invr,
 401                                                     const float*               E_lj)
 402 {
 403     /* potential switch constants */
 404     const float switch_V3 = nbparam->vdw_switch.c3;
 405     const float switch_V4 = nbparam->vdw_switch.c4;
 406     const float switch_V5 = nbparam->vdw_switch.c5;
 407     const float switch_F2 = nbparam->vdw_switch.c3;
 408     const float switch_F3 = nbparam->vdw_switch.c4;
 409     const float switch_F4 = nbparam->vdw_switch.c5;
 410
 411     const float r        = r2 * inv_r;
 412     const float r_switch = r - nbparam->rvdw_switch;
 413
 414     /* Unlike in the F+E kernel, conditional is faster here */
 415     if (r_switch > 0.0F)
 416     {
 417         const float sw = 1.0F
 418                          + (switch_V3 + (switch_V4 + switch_V5 * r_switch) * r_switch) * r_switch
 419                                    * r_switch * r_switch;
 420         const float dsw = (switch_F2 + (switch_F3 + switch_F4 * r_switch) * r_switch) * r_switch * r_switch;
 421
 422         *F_invr = (*F_invr) * sw - inv_r * (*E_lj) * dsw;
 423     }
 424 }
 425
 426 /*! Apply potential switch, force + energy version. */
 427 gmx_opencl_inline void calculate_potential_switch_F_E(const cl_nbparam_params_t* nbparam,
 428                                                       float                      inv_r,
 429                                                       float                      r2,
 430                                                       float*                     F_invr,
 431                                                       float*                     E_lj)
 432 {
 433     /* potential switch constants */
 434     const float switch_V3 = nbparam->vdw_switch.c3;
 435     const float switch_V4 = nbparam->vdw_switch.c4;
 436     const float switch_V5 = nbparam->vdw_switch.c5;
 437     const float switch_F2 = nbparam->vdw_switch.c3;
 438     const float switch_F3 = nbparam->vdw_switch.c4;
 439     const float switch_F4 = nbparam->vdw_switch.c5;
 440
 441     const float r        = r2 * inv_r;
 442     float       r_switch = r - nbparam->rvdw_switch;
 443     r_switch             = r_switch >= 0.0F ? r_switch : 0.0F;
 444
 445     /* Unlike in the F-only kernel, masking is faster here */
 446     const float sw =
 447             1.0F + (switch_V3 + (switch_V4 + switch_V5 * r_switch) * r_switch) * r_switch * r_switch * r_switch;
 448     const float dsw = (switch_F2 + (switch_F3 + switch_F4 * r_switch) * r_switch) * r_switch * r_switch;
 449
 450     *F_invr = (*F_invr) * sw - inv_r * (*E_lj) * dsw;
 451     *E_lj *= sw;
 452 }
 453
 454 /*! Calculate LJ-PME grid force contribution with
 455  *  geometric combination rule.
 456  */
 457 gmx_opencl_inline void calculate_lj_ewald_comb_geom_F(__constant const float* nbfp_comb_climg2d,
 458                                                       int                     typei,
 459                                                       int                     typej,
 460                                                       float                   r2,
 461                                                       float                   inv_r2,
 462                                                       float                   lje_coeff2,
 463                                                       float                   lje_coeff6_6,
 464                                                       float*                  F_invr)
 465 {
 466     const float c6grid = nbfp_comb_climg2d[2 * typei] * nbfp_comb_climg2d[2 * typej];
 467
 468     /* Recalculate inv_r6 without exclusion mask */
 469     const float inv_r6_nm = inv_r2 * inv_r2 * inv_r2;
 470     const float cr2       = lje_coeff2 * r2;
 471     const float expmcr2   = exp(-cr2);
 472     const float poly      = 1.0F + cr2 + HALF_F * cr2 * cr2;
 473
 474     /* Subtract the grid force from the total LJ force */
 475     *F_invr += c6grid * (inv_r6_nm - expmcr2 * (inv_r6_nm * poly + lje_coeff6_6)) * inv_r2;
 476 }
 477
 478 /*! Calculate LJ-PME grid force + energy contribution with
 479  *  geometric combination rule.
 480  */
 481 gmx_opencl_inline void calculate_lj_ewald_comb_geom_F_E(__constant const float* nbfp_comb_climg2d,
 482                                                         const cl_nbparam_params_t* nbparam,
 483                                                         int                        typei,
 484                                                         int                        typej,
 485                                                         float                      r2,
 486                                                         float                      inv_r2,
 487                                                         float                      lje_coeff2,
 488                                                         float                      lje_coeff6_6,
 489                                                         float                      int_bit,
 490                                                         float*                     F_invr,
 491                                                         float*                     E_lj)
 492 {
 493     const float c6grid = nbfp_comb_climg2d[2 * typei] * nbfp_comb_climg2d[2 * typej];
 494
 495     /* Recalculate inv_r6 without exclusion mask */
 496     const float inv_r6_nm = inv_r2 * inv_r2 * inv_r2;
 497     const float cr2       = lje_coeff2 * r2;
 498     const float expmcr2   = exp(-cr2);
 499     const float poly      = 1.0F + cr2 + HALF_F * cr2 * cr2;
 500
 501     /* Subtract the grid force from the total LJ force */
 502     *F_invr += c6grid * (inv_r6_nm - expmcr2 * (inv_r6_nm * poly + lje_coeff6_6)) * inv_r2;
 503
 504     /* Shift should be applied only to real LJ pairs */
 505     const float sh_mask = nbparam->sh_lj_ewald * int_bit;
 506     *E_lj += ONE_SIXTH_F * c6grid * (inv_r6_nm * (1.0F - expmcr2 * poly) + sh_mask);
 507 }
 508
 509 /*! Calculate LJ-PME grid force + energy contribution (if E_lj != NULL) with
 510  *  Lorentz-Berthelot combination rule.
 511  *  We use a single F+E kernel with conditional because the performance impact
 512  *  of this is pretty small and LB on the CPU is anyway very slow.
 513  */
 514 gmx_opencl_inline void calculate_lj_ewald_comb_LB_F_E(__constant const float*    nbfp_comb_climg2d,
 515                                                       const cl_nbparam_params_t* nbparam,
 516                                                       int                        typei,
 517                                                       int                        typej,
 518                                                       float                      r2,
 519                                                       float                      inv_r2,
 520                                                       float                      lje_coeff2,
 521                                                       float                      lje_coeff6_6,
 522                                                       float                      int_bit,
 523                                                       bool                       with_E_lj,
 524                                                       float*                     F_invr,
 525                                                       float*                     E_lj)
 526 {
 527     /* sigma and epsilon are scaled to give 6*C6 */
 528     const float sigma = nbfp_comb_climg2d[2 * typei] + nbfp_comb_climg2d[2 * typej];
 529
 530     const float epsilon = nbfp_comb_climg2d[2 * typei + 1] * nbfp_comb_climg2d[2 * typej + 1];
 531
 532     const float sigma2 = sigma * sigma;
 533     const float c6grid = epsilon * sigma2 * sigma2 * sigma2;
 534
 535     /* Recalculate inv_r6 without exclusion mask */
 536     const float inv_r6_nm = inv_r2 * inv_r2 * inv_r2;
 537     const float cr2       = lje_coeff2 * r2;
 538     const float expmcr2   = exp(-cr2);
 539     const float poly      = 1.0F + cr2 + HALF_F * cr2 * cr2;
 540
 541     /* Subtract the grid force from the total LJ force */
 542     *F_invr += c6grid * (inv_r6_nm - expmcr2 * (inv_r6_nm * poly + lje_coeff6_6)) * inv_r2;
 543
 544     if (with_E_lj)
 545     {
 546
 547         /* Shift should be applied only to real LJ pairs */
 548         const float sh_mask = nbparam->sh_lj_ewald * int_bit;
 549         *E_lj += ONE_SIXTH_F * c6grid * (inv_r6_nm * (1.0F - expmcr2 * poly) + sh_mask);
 550     }
 551 }
 552
 553 /*! Interpolate Ewald coulomb force using the table through the tex_nbfp texture.
 554  *  Original idea: from the OpenMM project
 555  */
 556 gmx_opencl_inline float interpolate_coulomb_force_r(__constant const float* coulomb_tab_climg2d,
 557                                                     float                   r,
 558                                                     float                   scale)
 559 {
 560     float normalized = scale * r;
 561     int   index      = (int)normalized;
 562     float fract2     = normalized - (float)index;
 563     float fract1     = 1.0F - fract2;
 564
 565     return fract1 * coulomb_tab_climg2d[index] + fract2 * coulomb_tab_climg2d[index + 1];
 566 }
 567
 568 /*! Calculate analytical Ewald correction term. */
 569 gmx_opencl_inline float pmecorrF(float z2)
 570 {
 571     const float FN6 = -1.7357322914161492954e-8F;
 572     const float FN5 = 1.4703624142580877519e-6F;
 573     const float FN4 = -0.000053401640219807709149F;
 574     const float FN3 = 0.0010054721316683106153F;
 575     const float FN2 = -0.019278317264888380590F;
 576     const float FN1 = 0.069670166153766424023F;
 577     const float FN0 = -0.75225204789749321333F;
 578
 579     const float FD4 = 0.0011193462567257629232F;
 580     const float FD3 = 0.014866955030185295499F;
 581     const float FD2 = 0.11583842382862377919F;
 582     const float FD1 = 0.50736591960530292870F;
 583     const float FD0 = 1.0F;
 584
 585     const float z4 = z2 * z2;
 586
 587     float polyFD0 = FD4 * z4 + FD2;
 588     float polyFD1 = FD3 * z4 + FD1;
 589     polyFD0       = polyFD0 * z4 + FD0;
 590     polyFD0       = polyFD1 * z2 + polyFD0;
 591
 592     polyFD0 = 1.0F / polyFD0;
 593
 594     float polyFN0 = FN6 * z4 + FN4;
 595     float polyFN1 = FN5 * z4 + FN3;
 596     polyFN0       = polyFN0 * z4 + FN2;
 597     polyFN1       = polyFN1 * z4 + FN1;
 598     polyFN0       = polyFN0 * z4 + FN0;
 599     polyFN0       = polyFN1 * z2 + polyFN0;
 600
 601     return polyFN0 * polyFD0;
 602 }
 603
 604 #    if REDUCE_SHUFFLE
 605 gmx_opencl_inline void
 606 reduce_force_j_shfl(float3 fin, __global float* fout, int gmx_unused tidxi, int gmx_unused tidxj, int aidx)
 607 {
 608     /* Only does reduction over 4 elements in cluster. Needs to be changed
 609      * for CL_SIZE>4. See CUDA code for required code */
 610     fin.x += intel_sub_group_shuffle_down(fin.x, fin.x, 1);
 611     fin.y += intel_sub_group_shuffle_up(fin.y, fin.y, 1);
 612     fin.z += intel_sub_group_shuffle_down(fin.z, fin.z, 1);
 613     if ((tidxi & 1) == 1)
 614     {
 615         fin.x = fin.y;
 616     }
 617     fin.x += intel_sub_group_shuffle_down(fin.x, fin.x, 2);
 618     fin.z += intel_sub_group_shuffle_up(fin.z, fin.z, 2);
 619     if (tidxi == 2)
 620     {
 621         fin.x = fin.z;
 622     }
 623     if (tidxi < 3)
 624     {
 625         atomicAdd_g_f(&fout[3 * aidx + tidxi], fin.x);
 626     }
 627 }
 628 #    endif
 629
 630 gmx_opencl_inline void
 631 reduce_force_j_generic(__local float* f_buf, float3 fcj_buf, __global float* fout, int tidxi, int tidxj, int aidx)
 632 {
 633     int tidx                      = tidxi + tidxj * CL_SIZE;
 634     f_buf[tidx]                   = fcj_buf.x;
 635     f_buf[FBUF_STRIDE + tidx]     = fcj_buf.y;
 636     f_buf[2 * FBUF_STRIDE + tidx] = fcj_buf.z;
 637
 638     /* Split the reduction between the first 3 column threads
 639        Threads with column id 0 will do the reduction for (float3).x components
 640        Threads with column id 1 will do the reduction for (float3).y components
 641        Threads with column id 2 will do the reduction for (float3).z components.
 642        The reduction is performed for each line tidxj of f_buf. */
 643     if (tidxi < 3)
 644     {
 645         float f = 0.0F;
 646         for (int j = tidxj * CL_SIZE; j < (tidxj + 1) * CL_SIZE; j++)
 647         {
 648             f += f_buf[FBUF_STRIDE * tidxi + j];
 649         }
 650
 651         atomicAdd_g_f(&fout[3 * aidx + tidxi], f);
 652     }
 653 }
 654
 655 /*! Final j-force reduction
 656  */
 657 gmx_opencl_inline void reduce_force_j(__local float gmx_unused* f_buf,
 658                                       float3                    fcj_buf,
 659                                       __global float*           fout,
 660                                       int                       tidxi,
 661                                       int                       tidxj,
 662                                       int                       aidx)
 663 {
 664 #    if REDUCE_SHUFFLE
 665     reduce_force_j_shfl(fcj_buf, fout, tidxi, tidxj, aidx);
 666 #    else
 667     reduce_force_j_generic(f_buf, fcj_buf, fout, tidxi, tidxj, aidx);
 668 #    endif
 669 }
 670
 671 #    if REDUCE_SHUFFLE
 672 gmx_opencl_inline void reduce_force_i_and_shift_shfl(float3*         fci_buf,
 673                                                      __global float* fout,
 674                                                      bool            bCalcFshift,
 675                                                      int             tidxi,
 676                                                      int             tidxj,
 677                                                      int             sci,
 678                                                      int             shift,
 679                                                      __global float* fshift)
 680 {
 681     /* Only does reduction over 4 elements in cluster (2 per warp). Needs to be changed
 682      * for CL_SIZE>4.*/
 683     float2 fshift_buf = 0;
 684     for (int ci_offset = 0; ci_offset < c_nbnxnGpuNumClusterPerSupercluster; ci_offset++)
 685     {
 686         int    aidx = (sci * c_nbnxnGpuNumClusterPerSupercluster + ci_offset) * CL_SIZE + tidxi;
 687         float3 fin  = fci_buf[ci_offset];
 688         fin.x += intel_sub_group_shuffle_down(fin.x, fin.x, CL_SIZE);
 689         fin.y += intel_sub_group_shuffle_up(fin.y, fin.y, CL_SIZE);
 690         fin.z += intel_sub_group_shuffle_down(fin.z, fin.z, CL_SIZE);
 691
 692         if (tidxj & 1)
 693         {
 694             fin.x = fin.y;
 695         }
 696         /* Threads 0,1 and 2,3 increment x,y for their warp */
 697         atomicAdd_g_f(&fout[3 * aidx + (tidxj & 1)], fin.x);
 698         if (bCalcFshift)
 699         {
 700             fshift_buf[0] += fin.x;
 701         }
 702         /* Threads 0 and 2 increment z for their warp */
 703         if ((tidxj & 1) == 0)
 704         {
 705             atomicAdd_g_f(&fout[3 * aidx + 2], fin.z);
 706             if (bCalcFshift)
 707             {
 708                 fshift_buf[1] += fin.z;
 709             }
 710         }
 711     }
 712     /* add up local shift forces into global mem */
 713     if (bCalcFshift)
 714     {
 715         // Threads 0,1 and 2,3 update x,y
 716         atomicAdd_g_f(&(fshift[3 * shift + (tidxj & 1)]), fshift_buf[0]);
 717         // Threads 0 and 2 update z
 718         if ((tidxj & 1) == 0)
 719         {
 720             atomicAdd_g_f(&(fshift[3 * shift + 2]), fshift_buf[1]);
 721         }
 722     }
 723 }
 724 #    endif
 725
 726 /*! Final i-force reduction; this implementation works only with power of two
 727  *  array sizes.
 728  */
 729 gmx_opencl_inline void reduce_force_i_and_shift_pow2(volatile __local float* f_buf,
 730                                                      float3*                 fci_buf,
 731                                                      __global float*         fout,
 732                                                      bool                    bCalcFshift,
 733                                                      int                     tidxi,
 734                                                      int                     tidxj,
 735                                                      int                     sci,
 736                                                      int                     shift,
 737                                                      __global float*         fshift)
 738 {
 739     float fshift_buf = 0;
 740     for (int ci_offset = 0; ci_offset < c_nbnxnGpuNumClusterPerSupercluster; ci_offset++)
 741     {
 742         int aidx = (sci * c_nbnxnGpuNumClusterPerSupercluster + ci_offset) * CL_SIZE + tidxi;
 743         int tidx = tidxi + tidxj * CL_SIZE;
 744         /* store i forces in shmem */
 745         f_buf[tidx]                   = fci_buf[ci_offset].x;
 746         f_buf[FBUF_STRIDE + tidx]     = fci_buf[ci_offset].y;
 747         f_buf[2 * FBUF_STRIDE + tidx] = fci_buf[ci_offset].z;
 748         barrier(CLK_LOCAL_MEM_FENCE);
 749
 750         /* Reduce the initial CL_SIZE values for each i atom to half
 751          * every step by using CL_SIZE * i threads.
 752          * Can't just use i as loop variable because than nvcc refuses to unroll.
 753          */
 754         int i = CL_SIZE / 2;
 755         for (int j = CL_SIZE_LOG2 - 1; j > 0; j--)
 756         {
 757             if (tidxj < i)
 758             {
 759
 760                 f_buf[tidxj * CL_SIZE + tidxi] += f_buf[(tidxj + i) * CL_SIZE + tidxi];
 761                 f_buf[FBUF_STRIDE + tidxj * CL_SIZE + tidxi] +=
 762                         f_buf[FBUF_STRIDE + (tidxj + i) * CL_SIZE + tidxi];
 763                 f_buf[2 * FBUF_STRIDE + tidxj * CL_SIZE + tidxi] +=
 764                         f_buf[2 * FBUF_STRIDE + (tidxj + i) * CL_SIZE + tidxi];
 765             }
 766             i >>= 1;
 767         }
 768         /* needed because
 769          * a) for CL_SIZE<8: id 2 (doing z in next block) is in 2nd warp
 770          * b) for all CL_SIZE a barrier is needed before f_buf is reused by next reduce_force_i call
 771          * TODO: Test on Nvidia for performance difference between having the barrier here or after the atomicAdd
 772          */
 773         barrier(CLK_LOCAL_MEM_FENCE);
 774
 775         /* i == 1, last reduction step, writing to global mem */
 776         /* Split the reduction between the first 3 line threads
 777            Threads with line id 0 will do the reduction for (float3).x components
 778            Threads with line id 1 will do the reduction for (float3).y components
 779            Threads with line id 2 will do the reduction for (float3).z components. */
 780         if (tidxj < 3)
 781         {
 782             float f = f_buf[tidxj * FBUF_STRIDE + tidxi] + f_buf[tidxj * FBUF_STRIDE + i * CL_SIZE + tidxi];
 783
 784             atomicAdd_g_f(&fout[3 * aidx + tidxj], f);
 785
 786             if (bCalcFshift)
 787             {
 788                 fshift_buf += f;
 789             }
 790         }
 791     }
 792     /* add up local shift forces into global mem */
 793     if (bCalcFshift)
 794     {
 795         /* Only threads with tidxj < 3 will update fshift.
 796            The threads performing the update, must be the same as the threads
 797            storing the reduction result above.
 798          */
 799         if (tidxj < 3)
 800         {
 801             atomicAdd_g_f(&(fshift[3 * shift + tidxj]), fshift_buf);
 802         }
 803     }
 804 }
 805
 806 /*! Final i-force reduction
 807  */
 808 gmx_opencl_inline void reduce_force_i_and_shift(__local float gmx_unused* f_buf,
 809                                                 float3*                   fci_buf,
 810                                                 __global float*           f,
 811                                                 bool                      bCalcFshift,
 812                                                 int                       tidxi,
 813                                                 int                       tidxj,
 814                                                 int                       sci,
 815                                                 int                       shift,
 816                                                 __global float*           fshift)
 817 {
 818 #    if REDUCE_SHUFFLE
 819     reduce_force_i_and_shift_shfl(fci_buf, f, bCalcFshift, tidxi, tidxj, sci, shift, fshift);
 820 #    else
 821     reduce_force_i_and_shift_pow2(f_buf, fci_buf, f, bCalcFshift, tidxi, tidxj, sci, shift, fshift);
 822 #    endif
 823 }
 824
 825
 826 #    if REDUCE_SHUFFLE
 827 gmx_opencl_inline void reduce_energy_shfl(float                    E_lj,
 828                                           float                    E_el,
 829                                           volatile __global float* e_lj,
 830                                           volatile __global float* e_el,
 831                                           int                      tidx)
 832 {
 833     E_lj = sub_group_reduce_add(E_lj);
 834     E_el = sub_group_reduce_add(E_el);
 835     /* Should be get_sub_group_local_id()==0. Doesn't work with Intel Classic driver.
 836      * To make it possible to use REDUCE_SHUFFLE with single subgroup per i-j pair
 837      * (e.g. subgroup size 16 with CL_SIZE 4), either this "if" needs to be changed or
 838      * the definition of WARP_SIZE (currently CL_SIZE*CL_SIZE/2) needs to be changed
 839      * (by supporting c_nbnxnGpuClusterpairSplit=1). */
 840     if (tidx == 0 || tidx == WARP_SIZE)
 841     {
 842         atomicAdd_g_f(e_lj, E_lj);
 843         atomicAdd_g_f(e_el, E_el);
 844     }
 845 }
 846 #    endif
 847
 848 /*! Energy reduction; this implementation works only with power of two
 849  *  array sizes.
 850  */
 851 gmx_opencl_inline void reduce_energy_pow2(volatile __local float*  buf,
 852                                           volatile __global float* e_lj,
 853                                           volatile __global float* e_el,
 854                                           int                      tidx)
 855 {
 856     int i = WARP_SIZE / 2;
 857
 858     /* Can't just use i as loop variable because than nvcc refuses to unroll. */
 859     for (int j = WARP_SIZE_LOG2 - 1; j > 0; j--)
 860     {
 861         if (tidx < i)
 862         {
 863             buf[tidx] += buf[tidx + i];
 864             buf[FBUF_STRIDE + tidx] += buf[FBUF_STRIDE + tidx + i];
 865         }
 866         i >>= 1;
 867     }
 868
 869     /* last reduction step, writing to global mem */
 870     if (tidx == 0)
 871     {
 872         float e1 = buf[tidx] + buf[tidx + i];
 873         float e2 = buf[FBUF_STRIDE + tidx] + buf[FBUF_STRIDE + tidx + i];
 874
 875         atomicAdd_g_f(e_lj, e1);
 876         atomicAdd_g_f(e_el, e2);
 877     }
 878 }
 879
 880 gmx_opencl_inline void reduce_energy(volatile __local float gmx_unused* buf,
 881                                      float                              E_lj,
 882                                      float                              E_el,
 883                                      volatile __global float*           e_lj,
 884                                      volatile __global float*           e_el,
 885                                      int                                tidx)
 886 {
 887 #    if REDUCE_SHUFFLE
 888     reduce_energy_shfl(E_lj, E_el, e_lj, e_el, tidx);
 889 #    else
 890     /* flush the energies to shmem and reduce them */
 891     buf[tidx] = E_lj;
 892     buf[FBUF_STRIDE + tidx] = E_el;
 893     reduce_energy_pow2(buf + (tidx & WARP_SIZE), e_lj, e_el, tidx & ~WARP_SIZE);
 894 #    endif
 895 }
 896
 897 gmx_opencl_inline bool gmx_sub_group_any_localmem(volatile __local int* warp_any, int widx, bool pred)
 898 {
 899     if (pred)
 900     {
 901         warp_any[widx] = 1;
 902     }
 903
 904     bool ret = warp_any[widx];
 905
 906     warp_any[widx] = 0;
 907
 908     return ret;
 909 }
 910
 911 //! Returns a true if predicate is true for any work item in warp
 912 gmx_opencl_inline bool gmx_sub_group_any(volatile __local int gmx_unused* warp_any, int gmx_unused widx, bool pred)
 913 {
 914 #    if USE_SUBGROUP_ANY
 915     return sub_group_any(pred);
 916 #    else
 917     return gmx_sub_group_any_localmem(warp_any, widx, pred);
 918 #    endif
 919 }
 920
 921 #endif /* NBNXN_OPENCL_KERNEL_UTILS_CLH */