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,2017,2018,2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \file
  36  *  \brief Define CUDA implementation of nbnxn_gpu.h
  37  *
  38  *  \author Szilard Pall <pall.szilard@gmail.com>
  39  */
  40 #include "gmxpre.h"
  41
  42 #include "config.h"
  43
  44 #include <assert.h>
  45 #include <stdlib.h>
  46
  47 #include "gromacs/nbnxm/nbnxm_gpu.h"
  48
  49 #if defined(_MSVC)
  50 #include <limits>
  51 #endif
  52
  53
  54 #include "nbnxm_cuda.h"
  55
  56 #include "gromacs/gpu_utils/cudautils.cuh"
  57 #include "gromacs/mdlib/force_flags.h"
  58 #include "gromacs/nbnxm/gpu_common.h"
  59 #include "gromacs/nbnxm/gpu_common_utils.h"
  60 #include "gromacs/nbnxm/gpu_data_mgmt.h"
  61 #include "gromacs/nbnxm/nbnxm.h"
  62 #include "gromacs/nbnxm/pairlist.h"
  63 #include "gromacs/timing/gpu_timing.h"
  64 #include "gromacs/utility/cstringutil.h"
  65 #include "gromacs/utility/gmxassert.h"
  66
  67 #include "nbnxm_cuda_types.h"
  68
  69 /***** The kernel declarations/definitions come here *****/
  70
  71 /* Top-level kernel declaration generation: will generate through multiple
  72  * inclusion the following flavors for all kernel declarations:
  73  * - force-only output;
  74  * - force and energy output;
  75  * - force-only with pair list pruning;
  76  * - force and energy output with pair list pruning.
  77  */
  78 #define FUNCTION_DECLARATION_ONLY
  79 /** Force only **/
  80 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  81 /** Force & energy **/
  82 #define CALC_ENERGIES
  83 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  84 #undef CALC_ENERGIES
  85
  86 /*** Pair-list pruning kernels ***/
  87 /** Force only **/
  88 #define PRUNE_NBL
  89 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  90 /** Force & energy **/
  91 #define CALC_ENERGIES
  92 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernels.cuh"
  93 #undef CALC_ENERGIES
  94 #undef PRUNE_NBL
  95
  96 /* Prune-only kernels */
  97 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_pruneonly.cuh"
  98 #undef FUNCTION_DECLARATION_ONLY
  99
 100 /* Now generate the function definitions if we are using a single compilation unit. */
 101 #if GMX_CUDA_NB_SINGLE_COMPILATION_UNIT
 102 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_F_noprune.cu"
 103 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_F_prune.cu"
 104 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_VF_noprune.cu"
 105 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_VF_prune.cu"
 106 #include "gromacs/nbnxm/cuda/nbnxm_cuda_kernel_pruneonly.cu"
 107 #endif /* GMX_CUDA_NB_SINGLE_COMPILATION_UNIT */
 108
 109
 110 /*! Nonbonded kernel function pointer type */
 111 typedef void (*nbnxn_cu_kfunc_ptr_t)(const cu_atomdata_t,
 112                                      const cu_nbparam_t,
 113                                      const cu_plist_t,
 114                                      bool);
 115
 116 /*********************************/
 117
 118 /*! Returns the number of blocks to be used for the nonbonded GPU kernel. */
 119 static inline int calc_nb_kernel_nblock(int nwork_units, const gmx_device_info_t *dinfo)
 120 {
 121     int max_grid_x_size;
 122
 123     assert(dinfo);
 124     /* CUDA does not accept grid dimension of 0 (which can happen e.g. with an
 125        empty domain) and that case should be handled before this point. */
 126     assert(nwork_units > 0);
 127
 128     max_grid_x_size = dinfo->prop.maxGridSize[0];
 129
 130     /* do we exceed the grid x dimension limit? */
 131     if (nwork_units > max_grid_x_size)
 132     {
 133         gmx_fatal(FARGS, "Watch out, the input system is too large to simulate!\n"
 134                   "The number of nonbonded work units (=number of super-clusters) exceeds the"
 135                   "maximum grid size in x dimension (%d > %d)!", nwork_units, max_grid_x_size);
 136     }
 137
 138     return nwork_units;
 139 }
 140
 141
 142 /* Constant arrays listing all kernel function pointers and enabling selection
 143    of a kernel in an elegant manner. */
 144
 145 /*! Pointers to the non-bonded kernels organized in 2-dim arrays by:
 146  *  electrostatics and VDW type.
 147  *
 148  *  Note that the row- and column-order of function pointers has to match the
 149  *  order of corresponding enumerated electrostatics and vdw types, resp.,
 150  *  defined in nbnxn_cuda_types.h.
 151  */
 152
 153 /*! Force-only kernel function pointers. */
 154 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_noprune_ptr[eelCuNR][evdwCuNR] =
 155 {
 156     { nbnxn_kernel_ElecCut_VdwLJ_F_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_F_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_F_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_F_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_F_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_F_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_F_cuda            },
 157     { nbnxn_kernel_ElecRF_VdwLJ_F_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_F_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_F_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_F_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_F_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_F_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_F_cuda             },
 158     { nbnxn_kernel_ElecEwQSTab_VdwLJ_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_F_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_F_cuda        },
 159     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_F_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_F_cuda },
 160     { nbnxn_kernel_ElecEw_VdwLJ_F_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_F_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_F_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_F_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_F_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_F_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_F_cuda             },
 161     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_F_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_F_cuda      }
 162 };
 163
 164 /*! Force + energy kernel function pointers. */
 165 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_noprune_ptr[eelCuNR][evdwCuNR] =
 166 {
 167     { nbnxn_kernel_ElecCut_VdwLJ_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_VF_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_VF_cuda            },
 168     { nbnxn_kernel_ElecRF_VdwLJ_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VF_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_VF_cuda             },
 169     { nbnxn_kernel_ElecEwQSTab_VdwLJ_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_VF_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_VF_cuda        },
 170     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_VF_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_VF_cuda },
 171     { nbnxn_kernel_ElecEw_VdwLJ_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VF_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_VF_cuda             },
 172     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VF_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_VF_cuda      }
 173 };
 174
 175 /*! Force + pruning kernel function pointers. */
 176 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_prune_ptr[eelCuNR][evdwCuNR] =
 177 {
 178     { nbnxn_kernel_ElecCut_VdwLJ_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_F_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_F_prune_cuda             },
 179     { nbnxn_kernel_ElecRF_VdwLJ_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_F_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_F_prune_cuda              },
 180     { nbnxn_kernel_ElecEwQSTab_VdwLJ_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_F_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_F_prune_cuda         },
 181     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_F_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_F_prune_cuda  },
 182     { nbnxn_kernel_ElecEw_VdwLJ_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_F_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_F_prune_cuda              },
 183     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_F_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_F_prune_cuda       }
 184 };
 185
 186 /*! Force + energy + pruning kernel function pointers. */
 187 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_prune_ptr[eelCuNR][evdwCuNR] =
 188 {
 189     { nbnxn_kernel_ElecCut_VdwLJ_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombGeom_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJCombLB_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJFsw_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJPsw_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombGeom_VF_prune_cuda,            nbnxn_kernel_ElecCut_VdwLJEwCombLB_VF_prune_cuda            },
 190     { nbnxn_kernel_ElecRF_VdwLJ_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJCombLB_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJFsw_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJPsw_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecRF_VdwLJEwCombLB_VF_prune_cuda             },
 191     { nbnxn_kernel_ElecEwQSTab_VdwLJ_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombGeom_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJCombLB_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJFsw_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJPsw_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombGeom_VF_prune_cuda,        nbnxn_kernel_ElecEwQSTab_VdwLJEwCombLB_VF_prune_cuda        },
 192     { nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJ_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombGeom_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJCombLB_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJFsw_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJPsw_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombGeom_VF_prune_cuda, nbnxn_kernel_ElecEwQSTabTwinCut_VdwLJEwCombLB_VF_prune_cuda },
 193     { nbnxn_kernel_ElecEw_VdwLJ_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJCombLB_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJFsw_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJPsw_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombGeom_VF_prune_cuda,             nbnxn_kernel_ElecEw_VdwLJEwCombLB_VF_prune_cuda             },
 194     { nbnxn_kernel_ElecEwTwinCut_VdwLJ_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombGeom_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJCombLB_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJFsw_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJPsw_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombGeom_VF_prune_cuda,      nbnxn_kernel_ElecEwTwinCut_VdwLJEwCombLB_VF_prune_cuda      }
 195 };
 196
 197 /*! Return a pointer to the kernel version to be executed at the current step. */
 198 static inline nbnxn_cu_kfunc_ptr_t select_nbnxn_kernel(int                                  eeltype,
 199                                                        int                                  evdwtype,
 200                                                        bool                                 bDoEne,
 201                                                        bool                                 bDoPrune,
 202                                                        const gmx_device_info_t gmx_unused  *devInfo)
 203 {
 204     nbnxn_cu_kfunc_ptr_t res;
 205
 206     GMX_ASSERT(eeltype < eelCuNR,
 207                "The electrostatics type requested is not implemented in the CUDA kernels.");
 208     GMX_ASSERT(evdwtype < evdwCuNR,
 209                "The VdW type requested is not implemented in the CUDA kernels.");
 210
 211     /* assert assumptions made by the kernels */
 212     GMX_ASSERT(c_nbnxnGpuClusterSize*c_nbnxnGpuClusterSize/c_nbnxnGpuClusterpairSplit == devInfo->prop.warpSize,
 213                "The CUDA kernels require the cluster_size_i*cluster_size_j/nbnxn_gpu_clusterpair_split to match the warp size of the architecture targeted.");
 214
 215     if (bDoEne)
 216     {
 217         if (bDoPrune)
 218         {
 219             res = nb_kfunc_ener_prune_ptr[eeltype][evdwtype];
 220         }
 221         else
 222         {
 223             res = nb_kfunc_ener_noprune_ptr[eeltype][evdwtype];
 224         }
 225     }
 226     else
 227     {
 228         if (bDoPrune)
 229         {
 230             res = nb_kfunc_noener_prune_ptr[eeltype][evdwtype];
 231         }
 232         else
 233         {
 234             res = nb_kfunc_noener_noprune_ptr[eeltype][evdwtype];
 235         }
 236     }
 237
 238     return res;
 239 }
 240
 241 /*! \brief Calculates the amount of shared memory required by the nonbonded kernel in use. */
 242 static inline int calc_shmem_required_nonbonded(const int num_threads_z, const gmx_device_info_t gmx_unused *dinfo, const cu_nbparam_t *nbp)
 243 {
 244     int shmem;
 245
 246     assert(dinfo);
 247
 248     /* size of shmem (force-buffers/xq/atom type preloading) */
 249     /* NOTE: with the default kernel on sm3.0 we need shmem only for pre-loading */
 250     /* i-atom x+q in shared memory */
 251     shmem  = c_numClPerSupercl * c_clSize * sizeof(float4);
 252     /* cj in shared memory, for each warp separately */
 253     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 254
 255     if (nbp->vdwtype == evdwCuCUTCOMBGEOM ||
 256         nbp->vdwtype == evdwCuCUTCOMBLB)
 257     {
 258         /* i-atom LJ combination parameters in shared memory */
 259         shmem += c_numClPerSupercl * c_clSize * sizeof(float2);
 260     }
 261     else
 262     {
 263         /* i-atom types in shared memory */
 264         shmem += c_numClPerSupercl * c_clSize * sizeof(int);
 265     }
 266
 267     return shmem;
 268 }
 269
 270 /*! \brief Launch asynchronously the xq buffer host to device copy. */
 271 void nbnxn_gpu_copy_xq_to_gpu(gmx_nbnxn_cuda_t       *nb,
 272                               const nbnxn_atomdata_t *nbatom,
 273                               int                     iloc,
 274                               bool                    haveOtherWork)
 275 {
 276     int                  adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 277
 278     cu_atomdata_t       *adat    = nb->atdat;
 279     cu_plist_t          *plist   = nb->plist[iloc];
 280     cu_timers_t         *t       = nb->timers;
 281     cudaStream_t         stream  = nb->stream[iloc];
 282
 283     bool                 bDoTime     = nb->bDoTime;
 284
 285     /* Don't launch the non-local H2D copy if there is no dependent
 286        work to do: neither non-local nor other (e.g. bonded) work
 287        to do that has as input the nbnxn coordaintes.
 288        Doing the same for the local kernel is more complicated, since the
 289        local part of the force array also depends on the non-local kernel.
 290        So to avoid complicating the code and to reduce the risk of bugs,
 291        we always call the local local x+q copy (and the rest of the local
 292        work in nbnxn_gpu_launch_kernel().
 293      */
 294     if (!haveOtherWork && canSkipWork(nb, iloc))
 295     {
 296         plist->haveFreshList = false;
 297
 298         return;
 299     }
 300
 301     /* calculate the atom data index range based on locality */
 302     if (LOCAL_I(iloc))
 303     {
 304         adat_begin  = 0;
 305         adat_len    = adat->natoms_local;
 306     }
 307     else
 308     {
 309         adat_begin  = adat->natoms_local;
 310         adat_len    = adat->natoms - adat->natoms_local;
 311     }
 312
 313     /* beginning of timed HtoD section */
 314     if (bDoTime)
 315     {
 316         t->nb_h2d[iloc].openTimingRegion(stream);
 317     }
 318
 319     /* HtoD x, q */
 320     cu_copy_H2D_async(adat->xq + adat_begin,
 321                       static_cast<const void *>(nbatom->x().data() + adat_begin * 4),
 322                       adat_len * sizeof(*adat->xq), stream);
 323
 324     if (bDoTime)
 325     {
 326         t->nb_h2d[iloc].closeTimingRegion(stream);
 327     }
 328
 329     /* When we get here all misc operations issued in the local stream as well as
 330        the local xq H2D are done,
 331        so we record that in the local stream and wait for it in the nonlocal one.
 332        This wait needs to precede any PP tasks, bonded or nonbonded, that may
 333        compute on interactions between local and nonlocal atoms.
 334      */
 335     if (nb->bUseTwoStreams)
 336     {
 337         if (iloc == eintLocal)
 338         {
 339             cudaError_t stat = cudaEventRecord(nb->misc_ops_and_local_H2D_done, stream);
 340             CU_RET_ERR(stat, "cudaEventRecord on misc_ops_and_local_H2D_done failed");
 341         }
 342         else
 343         {
 344             cudaError_t stat = cudaStreamWaitEvent(stream, nb->misc_ops_and_local_H2D_done, 0);
 345             CU_RET_ERR(stat, "cudaStreamWaitEvent on misc_ops_and_local_H2D_done failed");
 346         }
 347     }
 348 }
 349
 350 /*! As we execute nonbonded workload in separate streams, before launching
 351    the kernel we need to make sure that he following operations have completed:
 352    - atomdata allocation and related H2D transfers (every nstlist step);
 353    - pair list H2D transfer (every nstlist step);
 354    - shift vector H2D transfer (every nstlist step);
 355    - force (+shift force and energy) output clearing (every step).
 356
 357    These operations are issued in the local stream at the beginning of the step
 358    and therefore always complete before the local kernel launch. The non-local
 359    kernel is launched after the local on the same device/context hence it is
 360    inherently scheduled after the operations in the local stream (including the
 361    above "misc_ops") on pre-GK110 devices with single hardware queue, but on later
 362    devices with multiple hardware queues the dependency needs to be enforced.
 363    We use the misc_ops_and_local_H2D_done event to record the point where
 364    the local x+q H2D (and all preceding) tasks are complete and synchronize
 365    with this event in the non-local stream before launching the non-bonded kernel.
 366  */
 367 void nbnxn_gpu_launch_kernel(gmx_nbnxn_cuda_t       *nb,
 368                              int                     flags,
 369                              int                     iloc)
 370 {
 371     /* CUDA kernel launch-related stuff */
 372     int                  nblock;
 373     dim3                 dim_block, dim_grid;
 374     nbnxn_cu_kfunc_ptr_t nb_kernel = nullptr; /* fn pointer to the nonbonded kernel */
 375
 376     cu_atomdata_t       *adat    = nb->atdat;
 377     cu_nbparam_t        *nbp     = nb->nbparam;
 378     cu_plist_t          *plist   = nb->plist[iloc];
 379     cu_timers_t         *t       = nb->timers;
 380     cudaStream_t         stream  = nb->stream[iloc];
 381
 382     bool                 bCalcEner   = flags & GMX_FORCE_ENERGY;
 383     bool                 bCalcFshift = flags & GMX_FORCE_VIRIAL;
 384     bool                 bDoTime     = nb->bDoTime;
 385
 386     /* Don't launch the non-local kernel if there is no work to do.
 387        Doing the same for the local kernel is more complicated, since the
 388        local part of the force array also depends on the non-local kernel.
 389        So to avoid complicating the code and to reduce the risk of bugs,
 390        we always call the local kernel, and later (not in
 391        this function) the stream wait, local f copyback and the f buffer
 392        clearing. All these operations, except for the local interaction kernel,
 393        are needed for the non-local interactions. The skip of the local kernel
 394        call is taken care of later in this function. */
 395     if (canSkipWork(nb, iloc))
 396     {
 397         plist->haveFreshList = false;
 398
 399         return;
 400     }
 401
 402     if (nbp->useDynamicPruning && plist->haveFreshList)
 403     {
 404         /* Prunes for rlistOuter and rlistInner, sets plist->haveFreshList=false
 405            (TODO: ATM that's the way the timing accounting can distinguish between
 406            separate prune kernel and combined force+prune, maybe we need a better way?).
 407          */
 408         nbnxn_gpu_launch_kernel_pruneonly(nb, iloc, 1);
 409     }
 410
 411     if (plist->nsci == 0)
 412     {
 413         /* Don't launch an empty local kernel (not allowed with CUDA) */
 414         return;
 415     }
 416
 417     /* beginning of timed nonbonded calculation section */
 418     if (bDoTime)
 419     {
 420         t->nb_k[iloc].openTimingRegion(stream);
 421     }
 422
 423     /* get the pointer to the kernel flavor we need to use */
 424     nb_kernel = select_nbnxn_kernel(nbp->eeltype,
 425                                     nbp->vdwtype,
 426                                     bCalcEner,
 427                                     (plist->haveFreshList && !nb->timers->didPrune[iloc]),
 428                                     nb->dev_info);
 429
 430     /* Kernel launch config:
 431      * - The thread block dimensions match the size of i-clusters, j-clusters,
 432      *   and j-cluster concurrency, in x, y, and z, respectively.
 433      * - The 1D block-grid contains as many blocks as super-clusters.
 434      */
 435     int num_threads_z = 1;
 436     if (nb->dev_info->prop.major == 3 && nb->dev_info->prop.minor == 7)
 437     {
 438         num_threads_z = 2;
 439     }
 440     nblock    = calc_nb_kernel_nblock(plist->nsci, nb->dev_info);
 441
 442     KernelLaunchConfig config;
 443     config.blockSize[0]     = c_clSize;
 444     config.blockSize[1]     = c_clSize;
 445     config.blockSize[2]     = num_threads_z;
 446     config.gridSize[0]      = nblock;
 447     config.sharedMemorySize = calc_shmem_required_nonbonded(num_threads_z, nb->dev_info, nbp);
 448     config.stream           = stream;
 449
 450     if (debug)
 451     {
 452         fprintf(debug, "Non-bonded GPU launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 453                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 454                 "\tShMem: %zu\n",
 455                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 456                 config.gridSize[0], config.gridSize[1], plist->nsci*c_numClPerSupercl,
 457                 c_numClPerSupercl, plist->na_c,
 458                 config.sharedMemorySize);
 459     }
 460
 461     auto      *timingEvent = bDoTime ? t->nb_k[iloc].fetchNextEvent() : nullptr;
 462     const auto kernelArgs  = prepareGpuKernelArguments(nb_kernel, config, adat, nbp, plist, &bCalcFshift);
 463     launchGpuKernel(nb_kernel, config, timingEvent, "k_calc_nb", kernelArgs);
 464
 465     if (bDoTime)
 466     {
 467         t->nb_k[iloc].closeTimingRegion(stream);
 468     }
 469
 470     if (GMX_NATIVE_WINDOWS)
 471     {
 472         /* Windows: force flushing WDDM queue */
 473         cudaStreamQuery(stream);
 474     }
 475 }
 476
 477 /*! Calculates the amount of shared memory required by the CUDA kernel in use. */
 478 static inline int calc_shmem_required_prune(const int num_threads_z)
 479 {
 480     int shmem;
 481
 482     /* i-atom x in shared memory */
 483     shmem  = c_numClPerSupercl * c_clSize * sizeof(float4);
 484     /* cj in shared memory, for each warp separately */
 485     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 486
 487     return shmem;
 488 }
 489
 490 void nbnxn_gpu_launch_kernel_pruneonly(gmx_nbnxn_cuda_t       *nb,
 491                                        int                     iloc,
 492                                        int                     numParts)
 493 {
 494     cu_atomdata_t       *adat    = nb->atdat;
 495     cu_nbparam_t        *nbp     = nb->nbparam;
 496     cu_plist_t          *plist   = nb->plist[iloc];
 497     cu_timers_t         *t       = nb->timers;
 498     cudaStream_t         stream  = nb->stream[iloc];
 499
 500     bool                 bDoTime = nb->bDoTime;
 501
 502     if (plist->haveFreshList)
 503     {
 504         GMX_ASSERT(numParts == 1, "With first pruning we expect 1 part");
 505
 506         /* Set rollingPruningNumParts to signal that it is not set */
 507         plist->rollingPruningNumParts = 0;
 508         plist->rollingPruningPart     = 0;
 509     }
 510     else
 511     {
 512         if (plist->rollingPruningNumParts == 0)
 513         {
 514             plist->rollingPruningNumParts = numParts;
 515         }
 516         else
 517         {
 518             GMX_ASSERT(numParts == plist->rollingPruningNumParts, "It is not allowed to change numParts in between list generation steps");
 519         }
 520     }
 521
 522     /* Use a local variable for part and update in plist, so we can return here
 523      * without duplicating the part increment code.
 524      */
 525     int part = plist->rollingPruningPart;
 526
 527     plist->rollingPruningPart++;
 528     if (plist->rollingPruningPart >= plist->rollingPruningNumParts)
 529     {
 530         plist->rollingPruningPart = 0;
 531     }
 532
 533     /* Compute the number of list entries to prune in this pass */
 534     int numSciInPart = (plist->nsci - part)/numParts;
 535
 536     /* Don't launch the kernel if there is no work to do (not allowed with CUDA) */
 537     if (numSciInPart <= 0)
 538     {
 539         plist->haveFreshList = false;
 540
 541         return;
 542     }
 543
 544     GpuRegionTimer *timer = nullptr;
 545     if (bDoTime)
 546     {
 547         timer = &(plist->haveFreshList ? t->prune_k[iloc] : t->rollingPrune_k[iloc]);
 548     }
 549
 550     /* beginning of timed prune calculation section */
 551     if (bDoTime)
 552     {
 553         timer->openTimingRegion(stream);
 554     }
 555
 556     /* Kernel launch config:
 557      * - The thread block dimensions match the size of i-clusters, j-clusters,
 558      *   and j-cluster concurrency, in x, y, and z, respectively.
 559      * - The 1D block-grid contains as many blocks as super-clusters.
 560      */
 561     int                num_threads_z  = c_cudaPruneKernelJ4Concurrency;
 562     int                nblock         = calc_nb_kernel_nblock(numSciInPart, nb->dev_info);
 563     KernelLaunchConfig config;
 564     config.blockSize[0]     = c_clSize;
 565     config.blockSize[1]     = c_clSize;
 566     config.blockSize[2]     = num_threads_z;
 567     config.gridSize[0]      = nblock;
 568     config.sharedMemorySize = calc_shmem_required_prune(num_threads_z);
 569     config.stream           = stream;
 570
 571     if (debug)
 572     {
 573         fprintf(debug, "Pruning GPU kernel launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 574                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 575                 "\tShMem: %zu\n",
 576                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 577                 config.gridSize[0], config.gridSize[1], numSciInPart*c_numClPerSupercl,
 578                 c_numClPerSupercl, plist->na_c,
 579                 config.sharedMemorySize);
 580     }
 581
 582     auto          *timingEvent  = bDoTime ? timer->fetchNextEvent() : nullptr;
 583     constexpr char kernelName[] = "k_pruneonly";
 584     const auto    &kernel       = plist->haveFreshList ? nbnxn_kernel_prune_cuda<true> : nbnxn_kernel_prune_cuda<false>;
 585     const auto     kernelArgs   = prepareGpuKernelArguments(kernel, config, adat, nbp, plist, &numParts, &part);
 586     launchGpuKernel(kernel, config, timingEvent, kernelName, kernelArgs);
 587
 588     /* TODO: consider a more elegant way to track which kernel has been called
 589        (combined or separate 1st pass prune, rolling prune). */
 590     if (plist->haveFreshList)
 591     {
 592         plist->haveFreshList         = false;
 593         /* Mark that pruning has been done */
 594         nb->timers->didPrune[iloc] = true;
 595     }
 596     else
 597     {
 598         /* Mark that rolling pruning has been done */
 599         nb->timers->didRollingPrune[iloc] = true;
 600     }
 601
 602     if (bDoTime)
 603     {
 604         timer->closeTimingRegion(stream);
 605     }
 606
 607     if (GMX_NATIVE_WINDOWS)
 608     {
 609         /* Windows: force flushing WDDM queue */
 610         cudaStreamQuery(stream);
 611     }
 612 }
 613
 614 void nbnxn_gpu_launch_cpyback(gmx_nbnxn_cuda_t       *nb,
 615                               nbnxn_atomdata_t       *nbatom,
 616                               int                     flags,
 617                               int                     aloc,
 618                               bool                    haveOtherWork)
 619 {
 620     cudaError_t stat;
 621     int         adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 622
 623     /* determine interaction locality from atom locality */
 624     int              iloc = gpuAtomToInteractionLocality(aloc);
 625
 626     cu_atomdata_t   *adat    = nb->atdat;
 627     cu_timers_t     *t       = nb->timers;
 628     bool             bDoTime = nb->bDoTime;
 629     cudaStream_t     stream  = nb->stream[iloc];
 630
 631     bool             bCalcEner   = flags & GMX_FORCE_ENERGY;
 632     bool             bCalcFshift = flags & GMX_FORCE_VIRIAL;
 633
 634     /* don't launch non-local copy-back if there was no non-local work to do */
 635     if (!haveOtherWork && canSkipWork(nb, iloc))
 636     {
 637         return;
 638     }
 639
 640     getGpuAtomRange(adat, aloc, &adat_begin, &adat_len);
 641
 642     /* beginning of timed D2H section */
 643     if (bDoTime)
 644     {
 645         t->nb_d2h[iloc].openTimingRegion(stream);
 646     }
 647
 648     /* With DD the local D2H transfer can only start after the non-local
 649        kernel has finished. */
 650     if (iloc == eintLocal && nb->bUseTwoStreams)
 651     {
 652         stat = cudaStreamWaitEvent(stream, nb->nonlocal_done, 0);
 653         CU_RET_ERR(stat, "cudaStreamWaitEvent on nonlocal_done failed");
 654     }
 655
 656     /* DtoH f */
 657     cu_copy_D2H_async(nbatom->out[0].f.data() + adat_begin * 3, adat->f + adat_begin,
 658                       (adat_len)*sizeof(*adat->f), stream);
 659
 660     /* After the non-local D2H is launched the nonlocal_done event can be
 661        recorded which signals that the local D2H can proceed. This event is not
 662        placed after the non-local kernel because we want the non-local data
 663        back first. */
 664     if (iloc == eintNonlocal)
 665     {
 666         stat = cudaEventRecord(nb->nonlocal_done, stream);
 667         CU_RET_ERR(stat, "cudaEventRecord on nonlocal_done failed");
 668     }
 669
 670     /* only transfer energies in the local stream */
 671     if (LOCAL_I(iloc))
 672     {
 673         /* DtoH fshift */
 674         if (bCalcFshift)
 675         {
 676             cu_copy_D2H_async(nb->nbst.fshift, adat->fshift,
 677                               SHIFTS * sizeof(*nb->nbst.fshift), stream);
 678         }
 679
 680         /* DtoH energies */
 681         if (bCalcEner)
 682         {
 683             cu_copy_D2H_async(nb->nbst.e_lj, adat->e_lj,
 684                               sizeof(*nb->nbst.e_lj), stream);
 685             cu_copy_D2H_async(nb->nbst.e_el, adat->e_el,
 686                               sizeof(*nb->nbst.e_el), stream);
 687         }
 688     }
 689
 690     if (bDoTime)
 691     {
 692         t->nb_d2h[iloc].closeTimingRegion(stream);
 693     }
 694 }
 695
 696 void nbnxn_cuda_set_cacheconfig()
 697 {
 698     cudaError_t stat;
 699
 700     for (int i = 0; i < eelCuNR; i++)
 701     {
 702         for (int j = 0; j < evdwCuNR; j++)
 703         {
 704             /* Default kernel 32/32 kB Shared/L1 */
 705             cudaFuncSetCacheConfig(nb_kfunc_ener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 706             cudaFuncSetCacheConfig(nb_kfunc_ener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 707             cudaFuncSetCacheConfig(nb_kfunc_noener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 708             stat = cudaFuncSetCacheConfig(nb_kfunc_noener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 709             CU_RET_ERR(stat, "cudaFuncSetCacheConfig failed");
 710         }
 711     }
 712 }