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