src/gromacs/nbnxm/cuda/nbnxm_cuda.cu

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
   5  * Copyright (c) 2017,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 /*! \file
  37  *  \brief Define CUDA implementation of nbnxn_gpu.h
  38  *
  39  *  \author Szilard Pall <pall.szilard@gmail.com>
  40  */
  41 #include "gmxpre.h"
  42
  43 #include "config.h"
  44
  45 #include <assert.h>
  46 #include <stdlib.h>
  47
  48 #include "gromacs/nbnxm/nbnxm_gpu.h"
  49
  50 #if defined(_MSVC)
  51 #    include <limits>
  52 #endif
  53
  54
  55 #include "nbnxm_cuda.h"
  56
  57 #include "gromacs/gpu_utils/cudautils.cuh"
  58 #include "gromacs/gpu_utils/gpueventsynchronizer.cuh"
  59 #include "gromacs/gpu_utils/vectype_ops.cuh"
  60 #include "gromacs/mdtypes/simulation_workload.h"
  61 #include "gromacs/nbnxm/atomdata.h"
  62 #include "gromacs/nbnxm/gpu_common.h"
  63 #include "gromacs/nbnxm/gpu_common_utils.h"
  64 #include "gromacs/nbnxm/gpu_data_mgmt.h"
  65 #include "gromacs/nbnxm/grid.h"
  66 #include "gromacs/nbnxm/nbnxm.h"
  67 #include "gromacs/nbnxm/pairlist.h"
  68 #include "gromacs/timing/gpu_timing.h"
  69 #include "gromacs/utility/cstringutil.h"
  70 #include "gromacs/utility/gmxassert.h"
  71
  72 #include "nbnxm_buffer_ops_kernels.cuh"
  73 #include "nbnxm_cuda_types.h"
  74
  75 /***** The kernel declarations/definitions come here *****/
  76
  77 /* Top-level kernel declaration generation: will generate through multiple
  78  * inclusion the following flavors for all kernel declarations:
  79  * - force-only output;
  80  * - force and energy output;
  81  * - force-only with pair list pruning;
  82  * - force and energy output with pair list pruning.
  83  */
  84 #define FUNCTION_DECLARATION_ONLY
  85 /** Force only **/
  86 #include "nbnxm_cuda_kernels.cuh"
  87 /** Force & energy **/
  88 #define CALC_ENERGIES
  89 #include "nbnxm_cuda_kernels.cuh"
  90 #undef CALC_ENERGIES
  91
  92 /*** Pair-list pruning kernels ***/
  93 /** Force only **/
  94 #define PRUNE_NBL
  95 #include "nbnxm_cuda_kernels.cuh"
  96 /** Force & energy **/
  97 #define CALC_ENERGIES
  98 #include "nbnxm_cuda_kernels.cuh"
  99 #undef CALC_ENERGIES
 100 #undef PRUNE_NBL
 101
 102 /* Prune-only kernels */
 103 #include "nbnxm_cuda_kernel_pruneonly.cuh"
 104 #undef FUNCTION_DECLARATION_ONLY
 105
 106 /* Now generate the function definitions if we are using a single compilation unit. */
 107 #if GMX_CUDA_NB_SINGLE_COMPILATION_UNIT
 108 #    include "nbnxm_cuda_kernel_F_noprune.cu"
 109 #    include "nbnxm_cuda_kernel_F_prune.cu"
 110 #    include "nbnxm_cuda_kernel_VF_noprune.cu"
 111 #    include "nbnxm_cuda_kernel_VF_prune.cu"
 112 #    include "nbnxm_cuda_kernel_pruneonly.cu"
 113 #endif /* GMX_CUDA_NB_SINGLE_COMPILATION_UNIT */
 114
 115 namespace Nbnxm
 116 {
 117
 118 //! Number of CUDA threads in a block
 119 // TODO Optimize this through experimentation
 120 constexpr static int c_bufOpsThreadsPerBlock = 128;
 121
 122 /*! Nonbonded kernel function pointer type */
 123 typedef void (*nbnxn_cu_kfunc_ptr_t)(const cu_atomdata_t, const cu_nbparam_t, const cu_plist_t, bool);
 124
 125 /*********************************/
 126
 127 /*! Returns the number of blocks to be used for the nonbonded GPU kernel. */
 128 static inline int calc_nb_kernel_nblock(int nwork_units, const gmx_device_info_t* dinfo)
 129 {
 130     int max_grid_x_size;
 131
 132     assert(dinfo);
 133     /* CUDA does not accept grid dimension of 0 (which can happen e.g. with an
 134        empty domain) and that case should be handled before this point. */
 135     assert(nwork_units > 0);
 136
 137     max_grid_x_size = dinfo->prop.maxGridSize[0];
 138
 139     /* do we exceed the grid x dimension limit? */
 140     if (nwork_units > max_grid_x_size)
 141     {
 142         gmx_fatal(FARGS,
 143                   "Watch out, the input system is too large to simulate!\n"
 144                   "The number of nonbonded work units (=number of super-clusters) exceeds the"
 145                   "maximum grid size in x dimension (%d > %d)!",
 146                   nwork_units, max_grid_x_size);
 147     }
 148
 149     return nwork_units;
 150 }
 151
 152
 153 /* Constant arrays listing all kernel function pointers and enabling selection
 154    of a kernel in an elegant manner. */
 155
 156 /*! Pointers to the non-bonded kernels organized in 2-dim arrays by:
 157  *  electrostatics and VDW type.
 158  *
 159  *  Note that the row- and column-order of function pointers has to match the
 160  *  order of corresponding enumerated electrostatics and vdw types, resp.,
 161  *  defined in nbnxn_cuda_types.h.
 162  */
 163
 164 /*! Force-only kernel function pointers. */
 165 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_noprune_ptr[eelCuNR][evdwCuNR] = {
 166     { nbnxn_kernel_ElecCut_VdwLJ_F_cuda, nbnxn_kernel_ElecCut_VdwLJCombGeom_F_cuda,
 167       nbnxn_kernel_ElecCut_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecCut_VdwLJFsw_F_cuda,
 168       nbnxn_kernel_ElecCut_VdwLJPsw_F_cuda, nbnxn_kernel_ElecCut_VdwLJEwCombGeom_F_cuda,
 169       nbnxn_kernel_ElecCut_VdwLJEwCombLB_F_cuda },
 170     { nbnxn_kernel_ElecRF_VdwLJ_F_cuda, nbnxn_kernel_ElecRF_VdwLJCombGeom_F_cuda,
 171       nbnxn_kernel_ElecRF_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecRF_VdwLJFsw_F_cuda,
 172       nbnxn_kernel_ElecRF_VdwLJPsw_F_cuda, nbnxn_kernel_ElecRF_VdwLJEwCombGeom_F_cuda,
 173       nbnxn_kernel_ElecRF_VdwLJEwCombLB_F_cuda },
 174     { nbnxn_kernel_ElecEwQSTab_VdwLJ_F_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_F_cuda,
 175       nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJFsw_F_cuda,
 176       nbnxn_kernel_ElecEwQSTab_VdwLJPsw_F_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_F_cuda,
 177       nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_F_cuda },
 178     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_F_cuda,
 179       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_F_cuda,
 180       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_F_cuda,
 181       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_F_cuda },
 182     { nbnxn_kernel_ElecEw_VdwLJ_F_cuda, nbnxn_kernel_ElecEw_VdwLJCombGeom_F_cuda,
 183       nbnxn_kernel_ElecEw_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecEw_VdwLJFsw_F_cuda,
 184       nbnxn_kernel_ElecEw_VdwLJPsw_F_cuda, nbnxn_kernel_ElecEw_VdwLJEwCombGeom_F_cuda,
 185       nbnxn_kernel_ElecEw_VdwLJEwCombLB_F_cuda },
 186     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_F_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_F_cuda,
 187       nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_F_cuda,
 188       nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_F_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_F_cuda,
 189       nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_F_cuda }
 190 };
 191
 192 /*! Force + energy kernel function pointers. */
 193 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_noprune_ptr[eelCuNR][evdwCuNR] = {
 194     { nbnxn_kernel_ElecCut_VdwLJ_VF_cuda, nbnxn_kernel_ElecCut_VdwLJCombGeom_VF_cuda,
 195       nbnxn_kernel_ElecCut_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecCut_VdwLJFsw_VF_cuda,
 196       nbnxn_kernel_ElecCut_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecCut_VdwLJEwCombGeom_VF_cuda,
 197       nbnxn_kernel_ElecCut_VdwLJEwCombLB_VF_cuda },
 198     { nbnxn_kernel_ElecRF_VdwLJ_VF_cuda, nbnxn_kernel_ElecRF_VdwLJCombGeom_VF_cuda,
 199       nbnxn_kernel_ElecRF_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecRF_VdwLJFsw_VF_cuda,
 200       nbnxn_kernel_ElecRF_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VF_cuda,
 201       nbnxn_kernel_ElecRF_VdwLJEwCombLB_VF_cuda },
 202     { nbnxn_kernel_ElecEwQSTab_VdwLJ_VF_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_VF_cuda,
 203       nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJFsw_VF_cuda,
 204       nbnxn_kernel_ElecEwQSTab_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_VF_cuda,
 205       nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_VF_cuda },
 206     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_VF_cuda,
 207       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_VF_cuda,
 208       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_VF_cuda,
 209       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_VF_cuda },
 210     { nbnxn_kernel_ElecEw_VdwLJ_VF_cuda, nbnxn_kernel_ElecEw_VdwLJCombGeom_VF_cuda,
 211       nbnxn_kernel_ElecEw_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecEw_VdwLJFsw_VF_cuda,
 212       nbnxn_kernel_ElecEw_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VF_cuda,
 213       nbnxn_kernel_ElecEw_VdwLJEwCombLB_VF_cuda },
 214     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_VF_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VF_cuda,
 215       nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_VF_cuda,
 216       nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VF_cuda,
 217       nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_VF_cuda }
 218 };
 219
 220 /*! Force + pruning kernel function pointers. */
 221 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_prune_ptr[eelCuNR][evdwCuNR] = {
 222     { nbnxn_kernel_ElecCut_VdwLJ_F_prune_cuda, nbnxn_kernel_ElecCut_VdwLJCombGeom_F_prune_cuda,
 223       nbnxn_kernel_ElecCut_VdwLJCombLB_F_prune_cuda, nbnxn_kernel_ElecCut_VdwLJFsw_F_prune_cuda,
 224       nbnxn_kernel_ElecCut_VdwLJPsw_F_prune_cuda, nbnxn_kernel_ElecCut_VdwLJEwCombGeom_F_prune_cuda,
 225       nbnxn_kernel_ElecCut_VdwLJEwCombLB_F_prune_cuda },
 226     { nbnxn_kernel_ElecRF_VdwLJ_F_prune_cuda, nbnxn_kernel_ElecRF_VdwLJCombGeom_F_prune_cuda,
 227       nbnxn_kernel_ElecRF_VdwLJCombLB_F_prune_cuda, nbnxn_kernel_ElecRF_VdwLJFsw_F_prune_cuda,
 228       nbnxn_kernel_ElecRF_VdwLJPsw_F_prune_cuda, nbnxn_kernel_ElecRF_VdwLJEwCombGeom_F_prune_cuda,
 229       nbnxn_kernel_ElecRF_VdwLJEwCombLB_F_prune_cuda },
 230     { nbnxn_kernel_ElecEwQSTab_VdwLJ_F_prune_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_F_prune_cuda,
 231       nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_F_prune_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJFsw_F_prune_cuda,
 232       nbnxn_kernel_ElecEwQSTab_VdwLJPsw_F_prune_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_F_prune_cuda,
 233       nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_F_prune_cuda },
 234     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_F_prune_cuda,
 235       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_F_prune_cuda,
 236       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_F_prune_cuda,
 237       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_F_prune_cuda,
 238       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_F_prune_cuda,
 239       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_F_prune_cuda },
 240     { nbnxn_kernel_ElecEw_VdwLJ_F_prune_cuda, nbnxn_kernel_ElecEw_VdwLJCombGeom_F_prune_cuda,
 241       nbnxn_kernel_ElecEw_VdwLJCombLB_F_prune_cuda, nbnxn_kernel_ElecEw_VdwLJFsw_F_prune_cuda,
 242       nbnxn_kernel_ElecEw_VdwLJPsw_F_prune_cuda, nbnxn_kernel_ElecEw_VdwLJEwCombGeom_F_prune_cuda,
 243       nbnxn_kernel_ElecEw_VdwLJEwCombLB_F_prune_cuda },
 244     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_F_prune_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_F_prune_cuda,
 245       nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_F_prune_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_F_prune_cuda,
 246       nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_F_prune_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_F_prune_cuda,
 247       nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_F_prune_cuda }
 248 };
 249
 250 /*! Force + energy + pruning kernel function pointers. */
 251 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_prune_ptr[eelCuNR][evdwCuNR] = {
 252     { nbnxn_kernel_ElecCut_VdwLJ_VF_prune_cuda, nbnxn_kernel_ElecCut_VdwLJCombGeom_VF_prune_cuda,
 253       nbnxn_kernel_ElecCut_VdwLJCombLB_VF_prune_cuda, nbnxn_kernel_ElecCut_VdwLJFsw_VF_prune_cuda,
 254       nbnxn_kernel_ElecCut_VdwLJPsw_VF_prune_cuda, nbnxn_kernel_ElecCut_VdwLJEwCombGeom_VF_prune_cuda,
 255       nbnxn_kernel_ElecCut_VdwLJEwCombLB_VF_prune_cuda },
 256     { nbnxn_kernel_ElecRF_VdwLJ_VF_prune_cuda, nbnxn_kernel_ElecRF_VdwLJCombGeom_VF_prune_cuda,
 257       nbnxn_kernel_ElecRF_VdwLJCombLB_VF_prune_cuda, nbnxn_kernel_ElecRF_VdwLJFsw_VF_prune_cuda,
 258       nbnxn_kernel_ElecRF_VdwLJPsw_VF_prune_cuda, nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VF_prune_cuda,
 259       nbnxn_kernel_ElecRF_VdwLJEwCombLB_VF_prune_cuda },
 260     { nbnxn_kernel_ElecEwQSTab_VdwLJ_VF_prune_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_VF_prune_cuda,
 261       nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_VF_prune_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJFsw_VF_prune_cuda,
 262       nbnxn_kernel_ElecEwQSTab_VdwLJPsw_VF_prune_cuda, nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_VF_prune_cuda,
 263       nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_VF_prune_cuda },
 264     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_VF_prune_cuda,
 265       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_VF_prune_cuda,
 266       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_VF_prune_cuda,
 267       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_VF_prune_cuda,
 268       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_VF_prune_cuda,
 269       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_VF_prune_cuda,
 270       nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_VF_prune_cuda },
 271     { nbnxn_kernel_ElecEw_VdwLJ_VF_prune_cuda, nbnxn_kernel_ElecEw_VdwLJCombGeom_VF_prune_cuda,
 272       nbnxn_kernel_ElecEw_VdwLJCombLB_VF_prune_cuda, nbnxn_kernel_ElecEw_VdwLJFsw_VF_prune_cuda,
 273       nbnxn_kernel_ElecEw_VdwLJPsw_VF_prune_cuda, nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VF_prune_cuda,
 274       nbnxn_kernel_ElecEw_VdwLJEwCombLB_VF_prune_cuda },
 275     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_VF_prune_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VF_prune_cuda,
 276       nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VF_prune_cuda,
 277       nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_VF_prune_cuda, nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_VF_prune_cuda,
 278       nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VF_prune_cuda,
 279       nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_VF_prune_cuda }
 280 };
 281
 282 /*! Return a pointer to the kernel version to be executed at the current step. */
 283 static inline nbnxn_cu_kfunc_ptr_t select_nbnxn_kernel(int                     eeltype,
 284                                                        int                     evdwtype,
 285                                                        bool                    bDoEne,
 286                                                        bool                    bDoPrune,
 287                                                        const gmx_device_info_t gmx_unused* devInfo)
 288 {
 289     nbnxn_cu_kfunc_ptr_t res;
 290
 291     GMX_ASSERT(eeltype < eelCuNR,
 292                "The electrostatics type requested is not implemented in the CUDA kernels.");
 293     GMX_ASSERT(evdwtype < evdwCuNR,
 294                "The VdW type requested is not implemented in the CUDA kernels.");
 295
 296     /* assert assumptions made by the kernels */
 297     GMX_ASSERT(c_nbnxnGpuClusterSize * c_nbnxnGpuClusterSize / c_nbnxnGpuClusterpairSplit
 298                        == devInfo->prop.warpSize,
 299                "The CUDA kernels require the "
 300                "cluster_size_i*cluster_size_j/nbnxn_gpu_clusterpair_split to match the warp size "
 301                "of the architecture targeted.");
 302
 303     if (bDoEne)
 304     {
 305         if (bDoPrune)
 306         {
 307             res = nb_kfunc_ener_prune_ptr[eeltype][evdwtype];
 308         }
 309         else
 310         {
 311             res = nb_kfunc_ener_noprune_ptr[eeltype][evdwtype];
 312         }
 313     }
 314     else
 315     {
 316         if (bDoPrune)
 317         {
 318             res = nb_kfunc_noener_prune_ptr[eeltype][evdwtype];
 319         }
 320         else
 321         {
 322             res = nb_kfunc_noener_noprune_ptr[eeltype][evdwtype];
 323         }
 324     }
 325
 326     return res;
 327 }
 328
 329 /*! \brief Calculates the amount of shared memory required by the nonbonded kernel in use. */
 330 static inline int calc_shmem_required_nonbonded(const int               num_threads_z,
 331                                                 const gmx_device_info_t gmx_unused* dinfo,
 332                                                 const cu_nbparam_t*                 nbp)
 333 {
 334     int shmem;
 335
 336     assert(dinfo);
 337
 338     /* size of shmem (force-buffers/xq/atom type preloading) */
 339     /* NOTE: with the default kernel on sm3.0 we need shmem only for pre-loading */
 340     /* i-atom x+q in shared memory */
 341     shmem = c_nbnxnGpuNumClusterPerSupercluster * c_clSize * sizeof(float4);
 342     /* cj in shared memory, for each warp separately */
 343     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 344
 345     if (nbp->vdwtype == evdwCuCUTCOMBGEOM || nbp->vdwtype == evdwCuCUTCOMBLB)
 346     {
 347         /* i-atom LJ combination parameters in shared memory */
 348         shmem += c_nbnxnGpuNumClusterPerSupercluster * c_clSize * sizeof(float2);
 349     }
 350     else
 351     {
 352         /* i-atom types in shared memory */
 353         shmem += c_nbnxnGpuNumClusterPerSupercluster * c_clSize * sizeof(int);
 354     }
 355
 356     return shmem;
 357 }
 358
 359 /*! \brief Sync the nonlocal stream with dependent tasks in the local queue.
 360  *
 361  *  As the point where the local stream tasks can be considered complete happens
 362  *  at the same call point where the nonlocal stream should be synced with the
 363  *  the local, this function records the event if called with the local stream as
 364  *  argument and inserts in the GPU stream a wait on the event on the nonlocal.
 365  */
 366 void nbnxnInsertNonlocalGpuDependency(const NbnxmGpu* nb, const InteractionLocality interactionLocality)
 367 {
 368     cudaStream_t stream = nb->stream[interactionLocality];
 369
 370     /* When we get here all misc operations issued in the local stream as well as
 371        the local xq H2D are done,
 372        so we record that in the local stream and wait for it in the nonlocal one.
 373        This wait needs to precede any PP tasks, bonded or nonbonded, that may
 374        compute on interactions between local and nonlocal atoms.
 375      */
 376     if (nb->bUseTwoStreams)
 377     {
 378         if (interactionLocality == InteractionLocality::Local)
 379         {
 380             cudaError_t stat = cudaEventRecord(nb->misc_ops_and_local_H2D_done, stream);
 381             CU_RET_ERR(stat, "cudaEventRecord on misc_ops_and_local_H2D_done failed");
 382         }
 383         else
 384         {
 385             cudaError_t stat = cudaStreamWaitEvent(stream, nb->misc_ops_and_local_H2D_done, 0);
 386             CU_RET_ERR(stat, "cudaStreamWaitEvent on misc_ops_and_local_H2D_done failed");
 387         }
 388     }
 389 }
 390
 391 /*! \brief Launch asynchronously the xq buffer host to device copy. */
 392 void gpu_copy_xq_to_gpu(NbnxmGpu* nb, const nbnxn_atomdata_t* nbatom, const AtomLocality atomLocality)
 393 {
 394     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 395
 396     GMX_ASSERT(atomLocality == AtomLocality::Local || atomLocality == AtomLocality::NonLocal,
 397                "Only local and non-local xq transfers are supported");
 398
 399     const InteractionLocality iloc = gpuAtomToInteractionLocality(atomLocality);
 400
 401     int adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 402
 403     cu_atomdata_t* adat   = nb->atdat;
 404     cu_plist_t*    plist  = nb->plist[iloc];
 405     cu_timers_t*   t      = nb->timers;
 406     cudaStream_t   stream = nb->stream[iloc];
 407
 408     bool bDoTime = nb->bDoTime;
 409
 410     /* Don't launch the non-local H2D copy if there is no dependent
 411        work to do: neither non-local nor other (e.g. bonded) work
 412        to do that has as input the nbnxn coordaintes.
 413        Doing the same for the local kernel is more complicated, since the
 414        local part of the force array also depends on the non-local kernel.
 415        So to avoid complicating the code and to reduce the risk of bugs,
 416        we always call the local local x+q copy (and the rest of the local
 417        work in nbnxn_gpu_launch_kernel().
 418      */
 419     if ((iloc == InteractionLocality::NonLocal) && !haveGpuShortRangeWork(*nb, iloc))
 420     {
 421         plist->haveFreshList = false;
 422
 423         return;
 424     }
 425
 426     /* calculate the atom data index range based on locality */
 427     if (atomLocality == AtomLocality::Local)
 428     {
 429         adat_begin = 0;
 430         adat_len   = adat->natoms_local;
 431     }
 432     else
 433     {
 434         adat_begin = adat->natoms_local;
 435         adat_len   = adat->natoms - adat->natoms_local;
 436     }
 437
 438     /* HtoD x, q */
 439     /* beginning of timed HtoD section */
 440     if (bDoTime)
 441     {
 442         t->xf[atomLocality].nb_h2d.openTimingRegion(stream);
 443     }
 444
 445     cu_copy_H2D_async(adat->xq + adat_begin,
 446                       static_cast<const void*>(nbatom->x().data() + adat_begin * 4),
 447                       adat_len * sizeof(*adat->xq), stream);
 448
 449     if (bDoTime)
 450     {
 451         t->xf[atomLocality].nb_h2d.closeTimingRegion(stream);
 452     }
 453
 454     /* When we get here all misc operations issued in the local stream as well as
 455        the local xq H2D are done,
 456        so we record that in the local stream and wait for it in the nonlocal one.
 457        This wait needs to precede any PP tasks, bonded or nonbonded, that may
 458        compute on interactions between local and nonlocal atoms.
 459      */
 460     nbnxnInsertNonlocalGpuDependency(nb, iloc);
 461 }
 462
 463 /*! As we execute nonbonded workload in separate streams, before launching
 464    the kernel we need to make sure that he following operations have completed:
 465    - atomdata allocation and related H2D transfers (every nstlist step);
 466    - pair list H2D transfer (every nstlist step);
 467    - shift vector H2D transfer (every nstlist step);
 468    - force (+shift force and energy) output clearing (every step).
 469
 470    These operations are issued in the local stream at the beginning of the step
 471    and therefore always complete before the local kernel launch. The non-local
 472    kernel is launched after the local on the same device/context hence it is
 473    inherently scheduled after the operations in the local stream (including the
 474    above "misc_ops") on pre-GK110 devices with single hardware queue, but on later
 475    devices with multiple hardware queues the dependency needs to be enforced.
 476    We use the misc_ops_and_local_H2D_done event to record the point where
 477    the local x+q H2D (and all preceding) tasks are complete and synchronize
 478    with this event in the non-local stream before launching the non-bonded kernel.
 479  */
 480 void gpu_launch_kernel(NbnxmGpu* nb, const gmx::StepWorkload& stepWork, const InteractionLocality iloc)
 481 {
 482     cu_atomdata_t* adat   = nb->atdat;
 483     cu_nbparam_t*  nbp    = nb->nbparam;
 484     cu_plist_t*    plist  = nb->plist[iloc];
 485     cu_timers_t*   t      = nb->timers;
 486     cudaStream_t   stream = nb->stream[iloc];
 487
 488     bool bDoTime = nb->bDoTime;
 489
 490     /* Don't launch the non-local kernel if there is no work to do.
 491        Doing the same for the local kernel is more complicated, since the
 492        local part of the force array also depends on the non-local kernel.
 493        So to avoid complicating the code and to reduce the risk of bugs,
 494        we always call the local kernel, and later (not in
 495        this function) the stream wait, local f copyback and the f buffer
 496        clearing. All these operations, except for the local interaction kernel,
 497        are needed for the non-local interactions. The skip of the local kernel
 498        call is taken care of later in this function. */
 499     if (canSkipNonbondedWork(*nb, iloc))
 500     {
 501         plist->haveFreshList = false;
 502
 503         return;
 504     }
 505
 506     if (nbp->useDynamicPruning && plist->haveFreshList)
 507     {
 508         /* Prunes for rlistOuter and rlistInner, sets plist->haveFreshList=false
 509            (TODO: ATM that's the way the timing accounting can distinguish between
 510            separate prune kernel and combined force+prune, maybe we need a better way?).
 511          */
 512         gpu_launch_kernel_pruneonly(nb, iloc, 1);
 513     }
 514
 515     if (plist->nsci == 0)
 516     {
 517         /* Don't launch an empty local kernel (not allowed with CUDA) */
 518         return;
 519     }
 520
 521     /* beginning of timed nonbonded calculation section */
 522     if (bDoTime)
 523     {
 524         t->interaction[iloc].nb_k.openTimingRegion(stream);
 525     }
 526
 527     /* Kernel launch config:
 528      * - The thread block dimensions match the size of i-clusters, j-clusters,
 529      *   and j-cluster concurrency, in x, y, and z, respectively.
 530      * - The 1D block-grid contains as many blocks as super-clusters.
 531      */
 532     int num_threads_z = 1;
 533     if (nb->dev_info->prop.major == 3 && nb->dev_info->prop.minor == 7)
 534     {
 535         num_threads_z = 2;
 536     }
 537     int nblock = calc_nb_kernel_nblock(plist->nsci, nb->dev_info);
 538
 539
 540     KernelLaunchConfig config;
 541     config.blockSize[0]     = c_clSize;
 542     config.blockSize[1]     = c_clSize;
 543     config.blockSize[2]     = num_threads_z;
 544     config.gridSize[0]      = nblock;
 545     config.sharedMemorySize = calc_shmem_required_nonbonded(num_threads_z, nb->dev_info, nbp);
 546     config.stream           = stream;
 547
 548     if (debug)
 549     {
 550         fprintf(debug,
 551                 "Non-bonded GPU launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 552                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 553                 "\tShMem: %zu\n",
 554                 config.blockSize[0], config.blockSize[1], config.blockSize[2], config.gridSize[0],
 555                 config.gridSize[1], plist->nsci * c_nbnxnGpuNumClusterPerSupercluster,
 556                 c_nbnxnGpuNumClusterPerSupercluster, plist->na_c, config.sharedMemorySize);
 557     }
 558
 559     auto*      timingEvent = bDoTime ? t->interaction[iloc].nb_k.fetchNextEvent() : nullptr;
 560     const auto kernel      = select_nbnxn_kernel(
 561             nbp->eeltype, nbp->vdwtype, stepWork.computeEnergy,
 562             (plist->haveFreshList && !nb->timers->interaction[iloc].didPrune), nb->dev_info);
 563     const auto kernelArgs =
 564             prepareGpuKernelArguments(kernel, config, adat, nbp, plist, &stepWork.computeVirial);
 565     launchGpuKernel(kernel, config, timingEvent, "k_calc_nb", kernelArgs);
 566
 567     if (bDoTime)
 568     {
 569         t->interaction[iloc].nb_k.closeTimingRegion(stream);
 570     }
 571
 572     if (GMX_NATIVE_WINDOWS)
 573     {
 574         /* Windows: force flushing WDDM queue */
 575         cudaStreamQuery(stream);
 576     }
 577 }
 578
 579 /*! Calculates the amount of shared memory required by the CUDA kernel in use. */
 580 static inline int calc_shmem_required_prune(const int num_threads_z)
 581 {
 582     int shmem;
 583
 584     /* i-atom x in shared memory */
 585     shmem = c_nbnxnGpuNumClusterPerSupercluster * c_clSize * sizeof(float4);
 586     /* cj in shared memory, for each warp separately */
 587     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 588
 589     return shmem;
 590 }
 591
 592 void gpu_launch_kernel_pruneonly(NbnxmGpu* nb, const InteractionLocality iloc, const int numParts)
 593 {
 594     cu_atomdata_t* adat   = nb->atdat;
 595     cu_nbparam_t*  nbp    = nb->nbparam;
 596     cu_plist_t*    plist  = nb->plist[iloc];
 597     cu_timers_t*   t      = nb->timers;
 598     cudaStream_t   stream = nb->stream[iloc];
 599
 600     bool bDoTime = nb->bDoTime;
 601
 602     if (plist->haveFreshList)
 603     {
 604         GMX_ASSERT(numParts == 1, "With first pruning we expect 1 part");
 605
 606         /* Set rollingPruningNumParts to signal that it is not set */
 607         plist->rollingPruningNumParts = 0;
 608         plist->rollingPruningPart     = 0;
 609     }
 610     else
 611     {
 612         if (plist->rollingPruningNumParts == 0)
 613         {
 614             plist->rollingPruningNumParts = numParts;
 615         }
 616         else
 617         {
 618             GMX_ASSERT(numParts == plist->rollingPruningNumParts,
 619                        "It is not allowed to change numParts in between list generation steps");
 620         }
 621     }
 622
 623     /* Use a local variable for part and update in plist, so we can return here
 624      * without duplicating the part increment code.
 625      */
 626     int part = plist->rollingPruningPart;
 627
 628     plist->rollingPruningPart++;
 629     if (plist->rollingPruningPart >= plist->rollingPruningNumParts)
 630     {
 631         plist->rollingPruningPart = 0;
 632     }
 633
 634     /* Compute the number of list entries to prune in this pass */
 635     int numSciInPart = (plist->nsci - part) / numParts;
 636
 637     /* Don't launch the kernel if there is no work to do (not allowed with CUDA) */
 638     if (numSciInPart <= 0)
 639     {
 640         plist->haveFreshList = false;
 641
 642         return;
 643     }
 644
 645     GpuRegionTimer* timer = nullptr;
 646     if (bDoTime)
 647     {
 648         timer = &(plist->haveFreshList ? t->interaction[iloc].prune_k : t->interaction[iloc].rollingPrune_k);
 649     }
 650
 651     /* beginning of timed prune calculation section */
 652     if (bDoTime)
 653     {
 654         timer->openTimingRegion(stream);
 655     }
 656
 657     /* Kernel launch config:
 658      * - The thread block dimensions match the size of i-clusters, j-clusters,
 659      *   and j-cluster concurrency, in x, y, and z, respectively.
 660      * - The 1D block-grid contains as many blocks as super-clusters.
 661      */
 662     int                num_threads_z = c_cudaPruneKernelJ4Concurrency;
 663     int                nblock        = calc_nb_kernel_nblock(numSciInPart, nb->dev_info);
 664     KernelLaunchConfig config;
 665     config.blockSize[0]     = c_clSize;
 666     config.blockSize[1]     = c_clSize;
 667     config.blockSize[2]     = num_threads_z;
 668     config.gridSize[0]      = nblock;
 669     config.sharedMemorySize = calc_shmem_required_prune(num_threads_z);
 670     config.stream           = stream;
 671
 672     if (debug)
 673     {
 674         fprintf(debug,
 675                 "Pruning GPU kernel launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 676                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 677                 "\tShMem: %zu\n",
 678                 config.blockSize[0], config.blockSize[1], config.blockSize[2], config.gridSize[0],
 679                 config.gridSize[1], numSciInPart * c_nbnxnGpuNumClusterPerSupercluster,
 680                 c_nbnxnGpuNumClusterPerSupercluster, plist->na_c, config.sharedMemorySize);
 681     }
 682
 683     auto*          timingEvent  = bDoTime ? timer->fetchNextEvent() : nullptr;
 684     constexpr char kernelName[] = "k_pruneonly";
 685     const auto     kernel =
 686             plist->haveFreshList ? nbnxn_kernel_prune_cuda<true> : nbnxn_kernel_prune_cuda<false>;
 687     const auto kernelArgs = prepareGpuKernelArguments(kernel, config, adat, nbp, plist, &numParts, &part);
 688     launchGpuKernel(kernel, config, timingEvent, kernelName, kernelArgs);
 689
 690     /* TODO: consider a more elegant way to track which kernel has been called
 691        (combined or separate 1st pass prune, rolling prune). */
 692     if (plist->haveFreshList)
 693     {
 694         plist->haveFreshList = false;
 695         /* Mark that pruning has been done */
 696         nb->timers->interaction[iloc].didPrune = true;
 697     }
 698     else
 699     {
 700         /* Mark that rolling pruning has been done */
 701         nb->timers->interaction[iloc].didRollingPrune = true;
 702     }
 703
 704     if (bDoTime)
 705     {
 706         timer->closeTimingRegion(stream);
 707     }
 708
 709     if (GMX_NATIVE_WINDOWS)
 710     {
 711         /* Windows: force flushing WDDM queue */
 712         cudaStreamQuery(stream);
 713     }
 714 }
 715
 716 void gpu_launch_cpyback(NbnxmGpu*                nb,
 717                         nbnxn_atomdata_t*        nbatom,
 718                         const gmx::StepWorkload& stepWork,
 719                         const AtomLocality       atomLocality)
 720 {
 721     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 722
 723     cudaError_t stat;
 724     int         adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 725
 726     /* determine interaction locality from atom locality */
 727     const InteractionLocality iloc = gpuAtomToInteractionLocality(atomLocality);
 728
 729     /* extract the data */
 730     cu_atomdata_t* adat    = nb->atdat;
 731     cu_timers_t*   t       = nb->timers;
 732     bool           bDoTime = nb->bDoTime;
 733     cudaStream_t   stream  = nb->stream[iloc];
 734
 735     /* don't launch non-local copy-back if there was no non-local work to do */
 736     if ((iloc == InteractionLocality::NonLocal) && !haveGpuShortRangeWork(*nb, iloc))
 737     {
 738         return;
 739     }
 740
 741     getGpuAtomRange(adat, atomLocality, &adat_begin, &adat_len);
 742
 743     /* beginning of timed D2H section */
 744     if (bDoTime)
 745     {
 746         t->xf[atomLocality].nb_d2h.openTimingRegion(stream);
 747     }
 748
 749     /* With DD the local D2H transfer can only start after the non-local
 750        kernel has finished. */
 751     if (iloc == InteractionLocality::Local && nb->bUseTwoStreams)
 752     {
 753         stat = cudaStreamWaitEvent(stream, nb->nonlocal_done, 0);
 754         CU_RET_ERR(stat, "cudaStreamWaitEvent on nonlocal_done failed");
 755     }
 756
 757     /* DtoH f
 758      * Skip if buffer ops / reduction is offloaded to the GPU.
 759      */
 760     if (!stepWork.useGpuFBufferOps)
 761     {
 762         cu_copy_D2H_async(nbatom->out[0].f.data() + adat_begin * 3, adat->f + adat_begin,
 763                           (adat_len) * sizeof(*adat->f), stream);
 764     }
 765
 766     /* After the non-local D2H is launched the nonlocal_done event can be
 767        recorded which signals that the local D2H can proceed. This event is not
 768        placed after the non-local kernel because we want the non-local data
 769        back first. */
 770     if (iloc == InteractionLocality::NonLocal)
 771     {
 772         stat = cudaEventRecord(nb->nonlocal_done, stream);
 773         CU_RET_ERR(stat, "cudaEventRecord on nonlocal_done failed");
 774     }
 775
 776     /* only transfer energies in the local stream */
 777     if (iloc == InteractionLocality::Local)
 778     {
 779         /* DtoH fshift when virial is needed */
 780         if (stepWork.computeVirial)
 781         {
 782             cu_copy_D2H_async(nb->nbst.fshift, adat->fshift, SHIFTS * sizeof(*nb->nbst.fshift), stream);
 783         }
 784
 785         /* DtoH energies */
 786         if (stepWork.computeEnergy)
 787         {
 788             cu_copy_D2H_async(nb->nbst.e_lj, adat->e_lj, sizeof(*nb->nbst.e_lj), stream);
 789             cu_copy_D2H_async(nb->nbst.e_el, adat->e_el, sizeof(*nb->nbst.e_el), stream);
 790         }
 791     }
 792
 793     if (bDoTime)
 794     {
 795         t->xf[atomLocality].nb_d2h.closeTimingRegion(stream);
 796     }
 797 }
 798
 799 void cuda_set_cacheconfig()
 800 {
 801     cudaError_t stat;
 802
 803     for (int i = 0; i < eelCuNR; i++)
 804     {
 805         for (int j = 0; j < evdwCuNR; j++)
 806         {
 807             /* Default kernel 32/32 kB Shared/L1 */
 808             cudaFuncSetCacheConfig(nb_kfunc_ener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 809             cudaFuncSetCacheConfig(nb_kfunc_ener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 810             cudaFuncSetCacheConfig(nb_kfunc_noener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 811             stat = cudaFuncSetCacheConfig(nb_kfunc_noener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 812             CU_RET_ERR(stat, "cudaFuncSetCacheConfig failed");
 813         }
 814     }
 815 }
 816
 817 /* X buffer operations on GPU: performs conversion from rvec to nb format. */
 818 void nbnxn_gpu_x_to_nbat_x(const Nbnxm::Grid&        grid,
 819                            bool                      setFillerCoords,
 820                            NbnxmGpu*                 nb,
 821                            DeviceBuffer<float>       d_x,
 822                            GpuEventSynchronizer*     xReadyOnDevice,
 823                            const Nbnxm::AtomLocality locality,
 824                            int                       gridId,
 825                            int                       numColumnsMax)
 826 {
 827     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 828
 829     cu_atomdata_t* adat = nb->atdat;
 830
 831     const int                  numColumns      = grid.numColumns();
 832     const int                  cellOffset      = grid.cellOffset();
 833     const int                  numAtomsPerCell = grid.numAtomsPerCell();
 834     Nbnxm::InteractionLocality interactionLoc  = gpuAtomToInteractionLocality(locality);
 835
 836     cudaStream_t stream = nb->stream[interactionLoc];
 837
 838     int numAtoms = grid.srcAtomEnd() - grid.srcAtomBegin();
 839     // avoid empty kernel launch, skip to inserting stream dependency
 840     if (numAtoms != 0)
 841     {
 842         // TODO: This will only work with CUDA
 843         GMX_ASSERT(d_x, "Need a valid device pointer");
 844
 845         // ensure that coordinates are ready on the device before launching the kernel
 846         GMX_ASSERT(xReadyOnDevice, "Need a valid GpuEventSynchronizer object");
 847         xReadyOnDevice->enqueueWaitEvent(stream);
 848
 849         KernelLaunchConfig config;
 850         config.blockSize[0] = c_bufOpsThreadsPerBlock;
 851         config.blockSize[1] = 1;
 852         config.blockSize[2] = 1;
 853         config.gridSize[0] = (grid.numCellsColumnMax() * numAtomsPerCell + c_bufOpsThreadsPerBlock - 1)
 854                              / c_bufOpsThreadsPerBlock;
 855         config.gridSize[1] = numColumns;
 856         config.gridSize[2] = 1;
 857         GMX_ASSERT(config.gridSize[0] > 0,
 858                    "Can not have empty grid, early return above avoids this");
 859         config.sharedMemorySize = 0;
 860         config.stream           = stream;
 861
 862         auto kernelFn = setFillerCoords ? nbnxn_gpu_x_to_nbat_x_kernel<true>
 863                                         : nbnxn_gpu_x_to_nbat_x_kernel<false>;
 864         float4*    d_xq          = adat->xq;
 865         const int* d_atomIndices = nb->atomIndices;
 866         const int* d_cxy_na      = &nb->cxy_na[numColumnsMax * gridId];
 867         const int* d_cxy_ind     = &nb->cxy_ind[numColumnsMax * gridId];
 868         const auto kernelArgs =
 869                 prepareGpuKernelArguments(kernelFn, config, &numColumns, &d_xq, &d_x, &d_atomIndices,
 870                                           &d_cxy_na, &d_cxy_ind, &cellOffset, &numAtomsPerCell);
 871         launchGpuKernel(kernelFn, config, nullptr, "XbufferOps", kernelArgs);
 872     }
 873
 874     // TODO: note that this is not necessary when there are no local atoms, that is:
 875     // (numAtoms == 0 && interactionLoc == InteractionLocality::Local)
 876     // but for now we avoid that optimization
 877     nbnxnInsertNonlocalGpuDependency(nb, interactionLoc);
 878 }
 879
 880 /* F buffer operations on GPU: performs force summations and conversion from nb to rvec format.
 881  *
 882  * NOTE: When the total force device buffer is reallocated and its size increases, it is cleared in
 883  *       Local stream. Hence, if accumulateForce is true, NonLocal stream should start accumulating
 884  *       forces only after Local stream already done so.
 885  */
 886 void nbnxn_gpu_add_nbat_f_to_f(const AtomLocality                         atomLocality,
 887                                DeviceBuffer<float>                        totalForcesDevice,
 888                                NbnxmGpu*                                  nb,
 889                                void*                                      pmeForcesDevice,
 890                                gmx::ArrayRef<GpuEventSynchronizer* const> dependencyList,
 891                                int                                        atomStart,
 892                                int                                        numAtoms,
 893                                bool                                       useGpuFPmeReduction,
 894                                bool                                       accumulateForce)
 895 {
 896     GMX_ASSERT(nb, "Need a valid nbnxn_gpu object");
 897     GMX_ASSERT(numAtoms != 0, "Cannot call function with no atoms");
 898     GMX_ASSERT(totalForcesDevice, "Need a valid totalForcesDevice pointer");
 899
 900     const InteractionLocality iLocality = gpuAtomToInteractionLocality(atomLocality);
 901     cudaStream_t              stream    = nb->stream[iLocality];
 902     cu_atomdata_t*            adat      = nb->atdat;
 903
 904     size_t gmx_used_in_debug numDependency = static_cast<size_t>((useGpuFPmeReduction == true))
 905                                              + static_cast<size_t>((accumulateForce == true));
 906     GMX_ASSERT(numDependency >= dependencyList.size(),
 907                "Mismatching number of dependencies and call signature");
 908
 909     // Enqueue wait on all dependencies passed
 910     for (auto const synchronizer : dependencyList)
 911     {
 912         synchronizer->enqueueWaitEvent(stream);
 913     }
 914
 915     /* launch kernel */
 916
 917     KernelLaunchConfig config;
 918     config.blockSize[0] = c_bufOpsThreadsPerBlock;
 919     config.blockSize[1] = 1;
 920     config.blockSize[2] = 1;
 921     config.gridSize[0]  = ((numAtoms + 1) + c_bufOpsThreadsPerBlock - 1) / c_bufOpsThreadsPerBlock;
 922     config.gridSize[1]  = 1;
 923     config.gridSize[2]  = 1;
 924     config.sharedMemorySize = 0;
 925     config.stream           = stream;
 926
 927     auto kernelFn = accumulateForce ? nbnxn_gpu_add_nbat_f_to_f_kernel<true, false>
 928                                     : nbnxn_gpu_add_nbat_f_to_f_kernel<false, false>;
 929
 930     if (useGpuFPmeReduction)
 931     {
 932         GMX_ASSERT(pmeForcesDevice, "Need a valid pmeForcesDevice pointer");
 933         kernelFn = accumulateForce ? nbnxn_gpu_add_nbat_f_to_f_kernel<true, true>
 934                                    : nbnxn_gpu_add_nbat_f_to_f_kernel<false, true>;
 935     }
 936
 937     const float3* d_fNB    = adat->f;
 938     const float3* d_fPme   = (float3*)pmeForcesDevice;
 939     float3*       d_fTotal = (float3*)totalForcesDevice;
 940     const int*    d_cell   = nb->cell;
 941
 942     const auto kernelArgs = prepareGpuKernelArguments(kernelFn, config, &d_fNB, &d_fPme, &d_fTotal,
 943                                                       &d_cell, &atomStart, &numAtoms);
 944
 945     launchGpuKernel(kernelFn, config, nullptr, "FbufferOps", kernelArgs);
 946
 947     if (atomLocality == AtomLocality::Local)
 948     {
 949         GMX_ASSERT(nb->localFReductionDone != nullptr,
 950                    "localFReductionDone has to be a valid pointer");
 951         nb->localFReductionDone->markEvent(stream);
 952     }
 953 }
 954
 955 } // namespace Nbnxm