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