src/gromacs/mdlib/nbnxn_cuda/nbnxn_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/mdlib/nbnxn_gpu.h"
  48
  49 #if defined(_MSVC)
  50 #include <limits>
  51 #endif
  52
  53
  54 #include "nbnxn_cuda.h"
  55
  56 #include "gromacs/gpu_utils/cudautils.cuh"
  57 #include "gromacs/mdlib/force_flags.h"
  58 #include "gromacs/mdlib/nb_verlet.h"
  59 #include "gromacs/mdlib/nbnxn_gpu_common.h"
  60 #include "gromacs/mdlib/nbnxn_gpu_common_utils.h"
  61 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
  62 #include "gromacs/mdlib/nbnxn_pairlist.h"
  63 #include "gromacs/timing/gpu_timing.h"
  64 #include "gromacs/utility/cstringutil.h"
  65 #include "gromacs/utility/gmxassert.h"
  66
  67 #include "nbnxn_cuda_types.h"
  68
  69
  70 /***** The kernel declarations/definitions come here *****/
  71
  72 /* Top-level kernel declaration generation: will generate through multiple
  73  * inclusion the following flavors for all kernel declarations:
  74  * - force-only output;
  75  * - force and energy output;
  76  * - force-only with pair list pruning;
  77  * - force and energy output with pair list pruning.
  78  */
  79 #define FUNCTION_DECLARATION_ONLY
  80 /** Force only **/
  81 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernels.cuh"
  82 /** Force & energy **/
  83 #define CALC_ENERGIES
  84 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernels.cuh"
  85 #undef CALC_ENERGIES
  86
  87 /*** Pair-list pruning kernels ***/
  88 /** Force only **/
  89 #define PRUNE_NBL
  90 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernels.cuh"
  91 /** Force & energy **/
  92 #define CALC_ENERGIES
  93 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernels.cuh"
  94 #undef CALC_ENERGIES
  95 #undef PRUNE_NBL
  96
  97 /* Prune-only kernels */
  98 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_pruneonly.cuh"
  99 #undef FUNCTION_DECLARATION_ONLY
 100
 101 /* Now generate the function definitions if we are using a single compilation unit. */
 102 #if GMX_CUDA_NB_SINGLE_COMPILATION_UNIT
 103 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_F_noprune.cu"
 104 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_F_prune.cu"
 105 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_VF_noprune.cu"
 106 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_VF_prune.cu"
 107 #include "gromacs/mdlib/nbnxn_cuda/nbnxn_cuda_kernel_pruneonly.cu"
 108 #endif /* GMX_CUDA_NB_SINGLE_COMPILATION_UNIT */
 109
 110
 111 /*! Nonbonded kernel function pointer type */
 112 typedef void (*nbnxn_cu_kfunc_ptr_t)(const cu_atomdata_t,
 113                                      const cu_nbparam_t,
 114                                      const cu_plist_t,
 115                                      bool);
 116
 117 /*********************************/
 118
 119 /*! Returns the number of blocks to be used for the nonbonded GPU kernel. */
 120 static inline int calc_nb_kernel_nblock(int nwork_units, const gmx_device_info_t *dinfo)
 121 {
 122     int max_grid_x_size;
 123
 124     assert(dinfo);
 125     /* CUDA does not accept grid dimension of 0 (which can happen e.g. with an
 126        empty domain) and that case should be handled before this point. */
 127     assert(nwork_units > 0);
 128
 129     max_grid_x_size = dinfo->prop.maxGridSize[0];
 130
 131     /* do we exceed the grid x dimension limit? */
 132     if (nwork_units > max_grid_x_size)
 133     {
 134         gmx_fatal(FARGS, "Watch out, the input system is too large to simulate!\n"
 135                   "The number of nonbonded work units (=number of super-clusters) exceeds the"
 136                   "maximum grid size in x dimension (%d > %d)!", nwork_units, max_grid_x_size);
 137     }
 138
 139     return nwork_units;
 140 }
 141
 142
 143 /* Constant arrays listing all kernel function pointers and enabling selection
 144    of a kernel in an elegant manner. */
 145
 146 /*! Pointers to the non-bonded kernels organized in 2-dim arrays by:
 147  *  electrostatics and VDW type.
 148  *
 149  *  Note that the row- and column-order of function pointers has to match the
 150  *  order of corresponding enumerated electrostatics and vdw types, resp.,
 151  *  defined in nbnxn_cuda_types.h.
 152  */
 153
 154 /*! Force-only kernel function pointers. */
 155 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_noprune_ptr[eelCuNR][evdwCuNR] =
 156 {
 157     { 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            },
 158     { 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             },
 159     { 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        },
 160     { 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 },
 161     { 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             },
 162     { 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      }
 163 };
 164
 165 /*! Force + energy kernel function pointers. */
 166 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_noprune_ptr[eelCuNR][evdwCuNR] =
 167 {
 168     { 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            },
 169     { 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             },
 170     { 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        },
 171     { 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 },
 172     { 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             },
 173     { 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      }
 174 };
 175
 176 /*! Force + pruning kernel function pointers. */
 177 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_noener_prune_ptr[eelCuNR][evdwCuNR] =
 178 {
 179     { 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             },
 180     { 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              },
 181     { 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         },
 182     { 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  },
 183     { 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              },
 184     { 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       }
 185 };
 186
 187 /*! Force + energy + pruning kernel function pointers. */
 188 static const nbnxn_cu_kfunc_ptr_t nb_kfunc_ener_prune_ptr[eelCuNR][evdwCuNR] =
 189 {
 190     { 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            },
 191     { 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             },
 192     { 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        },
 193     { 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 },
 194     { 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             },
 195     { 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      }
 196 };
 197
 198 /*! Return a pointer to the kernel version to be executed at the current step. */
 199 static inline nbnxn_cu_kfunc_ptr_t select_nbnxn_kernel(int                                  eeltype,
 200                                                        int                                  evdwtype,
 201                                                        bool                                 bDoEne,
 202                                                        bool                                 bDoPrune,
 203                                                        const gmx_device_info_t gmx_unused  *devInfo)
 204 {
 205     nbnxn_cu_kfunc_ptr_t res;
 206
 207     GMX_ASSERT(eeltype < eelCuNR,
 208                "The electrostatics type requested is not implemented in the CUDA kernels.");
 209     GMX_ASSERT(evdwtype < evdwCuNR,
 210                "The VdW type requested is not implemented in the CUDA kernels.");
 211
 212     /* assert assumptions made by the kernels */
 213     GMX_ASSERT(c_nbnxnGpuClusterSize*c_nbnxnGpuClusterSize/c_nbnxnGpuClusterpairSplit == devInfo->prop.warpSize,
 214                "The CUDA kernels require the cluster_size_i*cluster_size_j/nbnxn_gpu_clusterpair_split to match the warp size of the architecture targeted.");
 215
 216     if (bDoEne)
 217     {
 218         if (bDoPrune)
 219         {
 220             res = nb_kfunc_ener_prune_ptr[eeltype][evdwtype];
 221         }
 222         else
 223         {
 224             res = nb_kfunc_ener_noprune_ptr[eeltype][evdwtype];
 225         }
 226     }
 227     else
 228     {
 229         if (bDoPrune)
 230         {
 231             res = nb_kfunc_noener_prune_ptr[eeltype][evdwtype];
 232         }
 233         else
 234         {
 235             res = nb_kfunc_noener_noprune_ptr[eeltype][evdwtype];
 236         }
 237     }
 238
 239     return res;
 240 }
 241
 242 /*! \brief Calculates the amount of shared memory required by the nonbonded kernel in use. */
 243 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)
 244 {
 245     int shmem;
 246
 247     assert(dinfo);
 248
 249     /* size of shmem (force-buffers/xq/atom type preloading) */
 250     /* NOTE: with the default kernel on sm3.0 we need shmem only for pre-loading */
 251     /* i-atom x+q in shared memory */
 252     shmem  = c_numClPerSupercl * c_clSize * sizeof(float4);
 253     /* cj in shared memory, for each warp separately */
 254     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 255
 256     if (nbp->vdwtype == evdwCuCUTCOMBGEOM ||
 257         nbp->vdwtype == evdwCuCUTCOMBLB)
 258     {
 259         /* i-atom LJ combination parameters in shared memory */
 260         shmem += c_numClPerSupercl * c_clSize * sizeof(float2);
 261     }
 262     else
 263     {
 264         /* i-atom types in shared memory */
 265         shmem += c_numClPerSupercl * c_clSize * sizeof(int);
 266     }
 267
 268     return shmem;
 269 }
 270
 271 /*! \brief Launch asynchronously the xq buffer host to device copy. */
 272 void nbnxn_gpu_copy_xq_to_gpu(gmx_nbnxn_cuda_t       *nb,
 273                               const nbnxn_atomdata_t *nbatom,
 274                               int                     iloc,
 275                               bool                    haveOtherWork)
 276 {
 277     int                  adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 278
 279     cu_atomdata_t       *adat    = nb->atdat;
 280     cu_plist_t          *plist   = nb->plist[iloc];
 281     cu_timers_t         *t       = nb->timers;
 282     cudaStream_t         stream  = nb->stream[iloc];
 283
 284     bool                 bDoTime     = nb->bDoTime;
 285
 286     /* Don't launch the non-local H2D copy if there is no dependent
 287        work to do: neither non-local nor other (e.g. bonded) work
 288        to do that has as input the nbnxn coordaintes.
 289        Doing the same for the local kernel is more complicated, since the
 290        local part of the force array also depends on the non-local kernel.
 291        So to avoid complicating the code and to reduce the risk of bugs,
 292        we always call the local local x+q copy (and the rest of the local
 293        work in nbnxn_gpu_launch_kernel().
 294      */
 295     if (!haveOtherWork && canSkipWork(nb, iloc))
 296     {
 297         plist->haveFreshList = false;
 298
 299         return;
 300     }
 301
 302     /* calculate the atom data index range based on locality */
 303     if (LOCAL_I(iloc))
 304     {
 305         adat_begin  = 0;
 306         adat_len    = adat->natoms_local;
 307     }
 308     else
 309     {
 310         adat_begin  = adat->natoms_local;
 311         adat_len    = adat->natoms - adat->natoms_local;
 312     }
 313
 314     /* beginning of timed HtoD section */
 315     if (bDoTime)
 316     {
 317         t->nb_h2d[iloc].openTimingRegion(stream);
 318     }
 319
 320     /* HtoD x, q */
 321     cu_copy_H2D_async(adat->xq + adat_begin,
 322                       static_cast<const void *>(nbatom->x().data() + adat_begin * 4),
 323                       adat_len * sizeof(*adat->xq), stream);
 324
 325     if (bDoTime)
 326     {
 327         t->nb_h2d[iloc].closeTimingRegion(stream);
 328     }
 329
 330     /* When we get here all misc operations issued in the local stream as well as
 331        the local xq H2D are done,
 332        so we record that in the local stream and wait for it in the nonlocal one.
 333        This wait needs to precede any PP tasks, bonded or nonbonded, that may
 334        compute on interactions between local and nonlocal atoms.
 335      */
 336     if (nb->bUseTwoStreams)
 337     {
 338         if (iloc == eintLocal)
 339         {
 340             cudaError_t stat = cudaEventRecord(nb->misc_ops_and_local_H2D_done, stream);
 341             CU_RET_ERR(stat, "cudaEventRecord on misc_ops_and_local_H2D_done failed");
 342         }
 343         else
 344         {
 345             cudaError_t stat = cudaStreamWaitEvent(stream, nb->misc_ops_and_local_H2D_done, 0);
 346             CU_RET_ERR(stat, "cudaStreamWaitEvent on misc_ops_and_local_H2D_done failed");
 347         }
 348     }
 349 }
 350
 351 /*! As we execute nonbonded workload in separate streams, before launching
 352    the kernel we need to make sure that he following operations have completed:
 353    - atomdata allocation and related H2D transfers (every nstlist step);
 354    - pair list H2D transfer (every nstlist step);
 355    - shift vector H2D transfer (every nstlist step);
 356    - force (+shift force and energy) output clearing (every step).
 357
 358    These operations are issued in the local stream at the beginning of the step
 359    and therefore always complete before the local kernel launch. The non-local
 360    kernel is launched after the local on the same device/context hence it is
 361    inherently scheduled after the operations in the local stream (including the
 362    above "misc_ops") on pre-GK110 devices with single hardware queue, but on later
 363    devices with multiple hardware queues the dependency needs to be enforced.
 364    We use the misc_ops_and_local_H2D_done event to record the point where
 365    the local x+q H2D (and all preceding) tasks are complete and synchronize
 366    with this event in the non-local stream before launching the non-bonded kernel.
 367  */
 368 void nbnxn_gpu_launch_kernel(gmx_nbnxn_cuda_t       *nb,
 369                              int                     flags,
 370                              int                     iloc)
 371 {
 372     /* CUDA kernel launch-related stuff */
 373     int                  nblock;
 374     dim3                 dim_block, dim_grid;
 375     nbnxn_cu_kfunc_ptr_t nb_kernel = nullptr; /* fn pointer to the nonbonded kernel */
 376
 377     cu_atomdata_t       *adat    = nb->atdat;
 378     cu_nbparam_t        *nbp     = nb->nbparam;
 379     cu_plist_t          *plist   = nb->plist[iloc];
 380     cu_timers_t         *t       = nb->timers;
 381     cudaStream_t         stream  = nb->stream[iloc];
 382
 383     bool                 bCalcEner   = flags & GMX_FORCE_ENERGY;
 384     bool                 bCalcFshift = flags & GMX_FORCE_VIRIAL;
 385     bool                 bDoTime     = nb->bDoTime;
 386
 387     /* Don't launch the non-local kernel if there is no work to do.
 388        Doing the same for the local kernel is more complicated, since the
 389        local part of the force array also depends on the non-local kernel.
 390        So to avoid complicating the code and to reduce the risk of bugs,
 391        we always call the local kernel, and later (not in
 392        this function) the stream wait, local f copyback and the f buffer
 393        clearing. All these operations, except for the local interaction kernel,
 394        are needed for the non-local interactions. The skip of the local kernel
 395        call is taken care of later in this function. */
 396     if (canSkipWork(nb, iloc))
 397     {
 398         plist->haveFreshList = false;
 399
 400         return;
 401     }
 402
 403     if (nbp->useDynamicPruning && plist->haveFreshList)
 404     {
 405         /* Prunes for rlistOuter and rlistInner, sets plist->haveFreshList=false
 406            (TODO: ATM that's the way the timing accounting can distinguish between
 407            separate prune kernel and combined force+prune, maybe we need a better way?).
 408          */
 409         nbnxn_gpu_launch_kernel_pruneonly(nb, iloc, 1);
 410     }
 411
 412     if (plist->nsci == 0)
 413     {
 414         /* Don't launch an empty local kernel (not allowed with CUDA) */
 415         return;
 416     }
 417
 418     /* beginning of timed nonbonded calculation section */
 419     if (bDoTime)
 420     {
 421         t->nb_k[iloc].openTimingRegion(stream);
 422     }
 423
 424     /* get the pointer to the kernel flavor we need to use */
 425     nb_kernel = select_nbnxn_kernel(nbp->eeltype,
 426                                     nbp->vdwtype,
 427                                     bCalcEner,
 428                                     (plist->haveFreshList && !nb->timers->didPrune[iloc]),
 429                                     nb->dev_info);
 430
 431     /* Kernel launch config:
 432      * - The thread block dimensions match the size of i-clusters, j-clusters,
 433      *   and j-cluster concurrency, in x, y, and z, respectively.
 434      * - The 1D block-grid contains as many blocks as super-clusters.
 435      */
 436     int num_threads_z = 1;
 437     if (nb->dev_info->prop.major == 3 && nb->dev_info->prop.minor == 7)
 438     {
 439         num_threads_z = 2;
 440     }
 441     nblock    = calc_nb_kernel_nblock(plist->nsci, nb->dev_info);
 442
 443     KernelLaunchConfig config;
 444     config.blockSize[0]     = c_clSize;
 445     config.blockSize[1]     = c_clSize;
 446     config.blockSize[2]     = num_threads_z;
 447     config.gridSize[0]      = nblock;
 448     config.sharedMemorySize = calc_shmem_required_nonbonded(num_threads_z, nb->dev_info, nbp);
 449     config.stream           = stream;
 450
 451     if (debug)
 452     {
 453         fprintf(debug, "Non-bonded GPU launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 454                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 455                 "\tShMem: %zu\n",
 456                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 457                 config.gridSize[0], config.gridSize[1], plist->nsci*c_numClPerSupercl,
 458                 c_numClPerSupercl, plist->na_c,
 459                 config.sharedMemorySize);
 460     }
 461
 462     auto      *timingEvent = bDoTime ? t->nb_k[iloc].fetchNextEvent() : nullptr;
 463     const auto kernelArgs  = prepareGpuKernelArguments(nb_kernel, config, adat, nbp, plist, &bCalcFshift);
 464     launchGpuKernel(nb_kernel, config, timingEvent, "k_calc_nb", kernelArgs);
 465
 466     if (bDoTime)
 467     {
 468         t->nb_k[iloc].closeTimingRegion(stream);
 469     }
 470
 471     if (GMX_NATIVE_WINDOWS)
 472     {
 473         /* Windows: force flushing WDDM queue */
 474         cudaStreamQuery(stream);
 475     }
 476 }
 477
 478 /*! Calculates the amount of shared memory required by the CUDA kernel in use. */
 479 static inline int calc_shmem_required_prune(const int num_threads_z)
 480 {
 481     int shmem;
 482
 483     /* i-atom x in shared memory */
 484     shmem  = c_numClPerSupercl * c_clSize * sizeof(float4);
 485     /* cj in shared memory, for each warp separately */
 486     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 487
 488     return shmem;
 489 }
 490
 491 void nbnxn_gpu_launch_kernel_pruneonly(gmx_nbnxn_cuda_t       *nb,
 492                                        int                     iloc,
 493                                        int                     numParts)
 494 {
 495     cu_atomdata_t       *adat    = nb->atdat;
 496     cu_nbparam_t        *nbp     = nb->nbparam;
 497     cu_plist_t          *plist   = nb->plist[iloc];
 498     cu_timers_t         *t       = nb->timers;
 499     cudaStream_t         stream  = nb->stream[iloc];
 500
 501     bool                 bDoTime = nb->bDoTime;
 502
 503     if (plist->haveFreshList)
 504     {
 505         GMX_ASSERT(numParts == 1, "With first pruning we expect 1 part");
 506
 507         /* Set rollingPruningNumParts to signal that it is not set */
 508         plist->rollingPruningNumParts = 0;
 509         plist->rollingPruningPart     = 0;
 510     }
 511     else
 512     {
 513         if (plist->rollingPruningNumParts == 0)
 514         {
 515             plist->rollingPruningNumParts = numParts;
 516         }
 517         else
 518         {
 519             GMX_ASSERT(numParts == plist->rollingPruningNumParts, "It is not allowed to change numParts in between list generation steps");
 520         }
 521     }
 522
 523     /* Use a local variable for part and update in plist, so we can return here
 524      * without duplicating the part increment code.
 525      */
 526     int part = plist->rollingPruningPart;
 527
 528     plist->rollingPruningPart++;
 529     if (plist->rollingPruningPart >= plist->rollingPruningNumParts)
 530     {
 531         plist->rollingPruningPart = 0;
 532     }
 533
 534     /* Compute the number of list entries to prune in this pass */
 535     int numSciInPart = (plist->nsci - part)/numParts;
 536
 537     /* Don't launch the kernel if there is no work to do (not allowed with CUDA) */
 538     if (numSciInPart <= 0)
 539     {
 540         plist->haveFreshList = false;
 541
 542         return;
 543     }
 544
 545     GpuRegionTimer *timer = nullptr;
 546     if (bDoTime)
 547     {
 548         timer = &(plist->haveFreshList ? t->prune_k[iloc] : t->rollingPrune_k[iloc]);
 549     }
 550
 551     /* beginning of timed prune calculation section */
 552     if (bDoTime)
 553     {
 554         timer->openTimingRegion(stream);
 555     }
 556
 557     /* Kernel launch config:
 558      * - The thread block dimensions match the size of i-clusters, j-clusters,
 559      *   and j-cluster concurrency, in x, y, and z, respectively.
 560      * - The 1D block-grid contains as many blocks as super-clusters.
 561      */
 562     int                num_threads_z  = c_cudaPruneKernelJ4Concurrency;
 563     int                nblock         = calc_nb_kernel_nblock(numSciInPart, nb->dev_info);
 564     KernelLaunchConfig config;
 565     config.blockSize[0]     = c_clSize;
 566     config.blockSize[1]     = c_clSize;
 567     config.blockSize[2]     = num_threads_z;
 568     config.gridSize[0]      = nblock;
 569     config.sharedMemorySize = calc_shmem_required_prune(num_threads_z);
 570     config.stream           = stream;
 571
 572     if (debug)
 573     {
 574         fprintf(debug, "Pruning GPU kernel launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 575                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 576                 "\tShMem: %zu\n",
 577                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 578                 config.gridSize[0], config.gridSize[1], numSciInPart*c_numClPerSupercl,
 579                 c_numClPerSupercl, plist->na_c,
 580                 config.sharedMemorySize);
 581     }
 582
 583     auto          *timingEvent  = bDoTime ? timer->fetchNextEvent() : nullptr;
 584     constexpr char kernelName[] = "k_pruneonly";
 585     const auto    &kernel       = plist->haveFreshList ? nbnxn_kernel_prune_cuda<true> : nbnxn_kernel_prune_cuda<false>;
 586     const auto     kernelArgs   = prepareGpuKernelArguments(kernel, config, adat, nbp, plist, &numParts, &part);
 587     launchGpuKernel(kernel, config, timingEvent, kernelName, kernelArgs);
 588
 589     /* TODO: consider a more elegant way to track which kernel has been called
 590        (combined or separate 1st pass prune, rolling prune). */
 591     if (plist->haveFreshList)
 592     {
 593         plist->haveFreshList         = false;
 594         /* Mark that pruning has been done */
 595         nb->timers->didPrune[iloc] = true;
 596     }
 597     else
 598     {
 599         /* Mark that rolling pruning has been done */
 600         nb->timers->didRollingPrune[iloc] = true;
 601     }
 602
 603     if (bDoTime)
 604     {
 605         timer->closeTimingRegion(stream);
 606     }
 607
 608     if (GMX_NATIVE_WINDOWS)
 609     {
 610         /* Windows: force flushing WDDM queue */
 611         cudaStreamQuery(stream);
 612     }
 613 }
 614
 615 void nbnxn_gpu_launch_cpyback(gmx_nbnxn_cuda_t       *nb,
 616                               nbnxn_atomdata_t       *nbatom,
 617                               int                     flags,
 618                               int                     aloc,
 619                               bool                    haveOtherWork)
 620 {
 621     cudaError_t stat;
 622     int         adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 623
 624     /* determine interaction locality from atom locality */
 625     int              iloc = gpuAtomToInteractionLocality(aloc);
 626
 627     cu_atomdata_t   *adat    = nb->atdat;
 628     cu_timers_t     *t       = nb->timers;
 629     bool             bDoTime = nb->bDoTime;
 630     cudaStream_t     stream  = nb->stream[iloc];
 631
 632     bool             bCalcEner   = flags & GMX_FORCE_ENERGY;
 633     bool             bCalcFshift = flags & GMX_FORCE_VIRIAL;
 634
 635     /* don't launch non-local copy-back if there was no non-local work to do */
 636     if (!haveOtherWork && canSkipWork(nb, iloc))
 637     {
 638         return;
 639     }
 640
 641     getGpuAtomRange(adat, aloc, &adat_begin, &adat_len);
 642
 643     /* beginning of timed D2H section */
 644     if (bDoTime)
 645     {
 646         t->nb_d2h[iloc].openTimingRegion(stream);
 647     }
 648
 649     /* With DD the local D2H transfer can only start after the non-local
 650        kernel has finished. */
 651     if (iloc == eintLocal && nb->bUseTwoStreams)
 652     {
 653         stat = cudaStreamWaitEvent(stream, nb->nonlocal_done, 0);
 654         CU_RET_ERR(stat, "cudaStreamWaitEvent on nonlocal_done failed");
 655     }
 656
 657     /* DtoH f */
 658     cu_copy_D2H_async(nbatom->out[0].f.data() + adat_begin * 3, adat->f + adat_begin,
 659                       (adat_len)*sizeof(*adat->f), stream);
 660
 661     /* After the non-local D2H is launched the nonlocal_done event can be
 662        recorded which signals that the local D2H can proceed. This event is not
 663        placed after the non-local kernel because we want the non-local data
 664        back first. */
 665     if (iloc == eintNonlocal)
 666     {
 667         stat = cudaEventRecord(nb->nonlocal_done, stream);
 668         CU_RET_ERR(stat, "cudaEventRecord on nonlocal_done failed");
 669     }
 670
 671     /* only transfer energies in the local stream */
 672     if (LOCAL_I(iloc))
 673     {
 674         /* DtoH fshift */
 675         if (bCalcFshift)
 676         {
 677             cu_copy_D2H_async(nb->nbst.fshift, adat->fshift,
 678                               SHIFTS * sizeof(*nb->nbst.fshift), stream);
 679         }
 680
 681         /* DtoH energies */
 682         if (bCalcEner)
 683         {
 684             cu_copy_D2H_async(nb->nbst.e_lj, adat->e_lj,
 685                               sizeof(*nb->nbst.e_lj), stream);
 686             cu_copy_D2H_async(nb->nbst.e_el, adat->e_el,
 687                               sizeof(*nb->nbst.e_el), stream);
 688         }
 689     }
 690
 691     if (bDoTime)
 692     {
 693         t->nb_d2h[iloc].closeTimingRegion(stream);
 694     }
 695 }
 696
 697 void nbnxn_cuda_set_cacheconfig()
 698 {
 699     cudaError_t stat;
 700
 701     for (int i = 0; i < eelCuNR; i++)
 702     {
 703         for (int j = 0; j < evdwCuNR; j++)
 704         {
 705             /* Default kernel 32/32 kB Shared/L1 */
 706             cudaFuncSetCacheConfig(nb_kfunc_ener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 707             cudaFuncSetCacheConfig(nb_kfunc_ener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 708             cudaFuncSetCacheConfig(nb_kfunc_noener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 709             stat = cudaFuncSetCacheConfig(nb_kfunc_noener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 710             CU_RET_ERR(stat, "cudaFuncSetCacheConfig failed");
 711         }
 712     }
 713 }