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, 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     if (dinfo->prop.major >= 3)
 256     {
 257         if (nbp->vdwtype == evdwCuCUTCOMBGEOM ||
 258             nbp->vdwtype == evdwCuCUTCOMBLB)
 259         {
 260             /* i-atom LJ combination parameters in shared memory */
 261             shmem += c_numClPerSupercl * c_clSize * sizeof(float2);
 262         }
 263         else
 264         {
 265             /* i-atom types in shared memory */
 266             shmem += c_numClPerSupercl * c_clSize * sizeof(int);
 267         }
 268     }
 269     if (dinfo->prop.major < 3)
 270     {
 271         /* force reduction buffers in shared memory */
 272         shmem += c_clSize * c_clSize * 3 * sizeof(float);
 273     }
 274     return shmem;
 275 }
 276
 277 /*! As we execute nonbonded workload in separate streams, before launching
 278    the kernel we need to make sure that he following operations have completed:
 279    - atomdata allocation and related H2D transfers (every nstlist step);
 280    - pair list H2D transfer (every nstlist step);
 281    - shift vector H2D transfer (every nstlist step);
 282    - force (+shift force and energy) output clearing (every step).
 283
 284    These operations are issued in the local stream at the beginning of the step
 285    and therefore always complete before the local kernel launch. The non-local
 286    kernel is launched after the local on the same device/context hence it is
 287    inherently scheduled after the operations in the local stream (including the
 288    above "misc_ops") on pre-GK110 devices with single hardware queue, but on later
 289    devices with multiple hardware queues the dependency needs to be enforced.
 290    We use the misc_ops_and_local_H2D_done event to record the point where
 291    the local x+q H2D (and all preceding) tasks are complete and synchronize
 292    with this event in the non-local stream before launching the non-bonded kernel.
 293  */
 294 void nbnxn_gpu_launch_kernel(gmx_nbnxn_cuda_t       *nb,
 295                              const nbnxn_atomdata_t *nbatom,
 296                              int                     flags,
 297                              int                     iloc)
 298 {
 299     cudaError_t          stat;
 300     int                  adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 301     /* CUDA kernel launch-related stuff */
 302     int                  nblock;
 303     dim3                 dim_block, dim_grid;
 304     nbnxn_cu_kfunc_ptr_t nb_kernel = NULL; /* fn pointer to the nonbonded kernel */
 305
 306     cu_atomdata_t       *adat    = nb->atdat;
 307     cu_nbparam_t        *nbp     = nb->nbparam;
 308     cu_plist_t          *plist   = nb->plist[iloc];
 309     cu_timers_t         *t       = nb->timers;
 310     cudaStream_t         stream  = nb->stream[iloc];
 311
 312     bool                 bCalcEner   = flags & GMX_FORCE_ENERGY;
 313     bool                 bCalcFshift = flags & GMX_FORCE_VIRIAL;
 314     bool                 bDoTime     = nb->bDoTime;
 315
 316     /* Don't launch the non-local kernel if there is no work to do.
 317        Doing the same for the local kernel is more complicated, since the
 318        local part of the force array also depends on the non-local kernel.
 319        So to avoid complicating the code and to reduce the risk of bugs,
 320        we always call the local kernel, the local x+q copy and later (not in
 321        this function) the stream wait, local f copyback and the f buffer
 322        clearing. All these operations, except for the local interaction kernel,
 323        are needed for the non-local interactions. The skip of the local kernel
 324        call is taken care of later in this function. */
 325     if (canSkipWork(nb, iloc))
 326     {
 327         plist->haveFreshList = false;
 328
 329         return;
 330     }
 331
 332     /* calculate the atom data index range based on locality */
 333     if (LOCAL_I(iloc))
 334     {
 335         adat_begin  = 0;
 336         adat_len    = adat->natoms_local;
 337     }
 338     else
 339     {
 340         adat_begin  = adat->natoms_local;
 341         adat_len    = adat->natoms - adat->natoms_local;
 342     }
 343
 344     /* beginning of timed HtoD section */
 345     if (bDoTime)
 346     {
 347         t->nb_h2d[iloc].openTimingRegion(stream);
 348     }
 349
 350     /* HtoD x, q */
 351     cu_copy_H2D_async(adat->xq + adat_begin, nbatom->x + adat_begin * 4,
 352                       adat_len * sizeof(*adat->xq), stream);
 353
 354     if (bDoTime)
 355     {
 356         t->nb_h2d[iloc].closeTimingRegion(stream);
 357     }
 358
 359     /* When we get here all misc operations issues in the local stream as well as
 360        the local xq H2D are done,
 361        so we record that in the local stream and wait for it in the nonlocal one. */
 362     if (nb->bUseTwoStreams)
 363     {
 364         if (iloc == eintLocal)
 365         {
 366             stat = cudaEventRecord(nb->misc_ops_and_local_H2D_done, stream);
 367             CU_RET_ERR(stat, "cudaEventRecord on misc_ops_and_local_H2D_done failed");
 368         }
 369         else
 370         {
 371             stat = cudaStreamWaitEvent(stream, nb->misc_ops_and_local_H2D_done, 0);
 372             CU_RET_ERR(stat, "cudaStreamWaitEvent on misc_ops_and_local_H2D_done failed");
 373         }
 374     }
 375
 376     if (nbp->useDynamicPruning && plist->haveFreshList)
 377     {
 378         /* Prunes for rlistOuter and rlistInner, sets plist->haveFreshList=false
 379            (TODO: ATM that's the way the timing accounting can distinguish between
 380            separate prune kernel and combined force+prune, maybe we need a better way?).
 381          */
 382         nbnxn_gpu_launch_kernel_pruneonly(nb, iloc, 1);
 383     }
 384
 385     if (plist->nsci == 0)
 386     {
 387         /* Don't launch an empty local kernel (not allowed with CUDA) */
 388         return;
 389     }
 390
 391     /* beginning of timed nonbonded calculation section */
 392     if (bDoTime)
 393     {
 394         t->nb_k[iloc].openTimingRegion(stream);
 395     }
 396
 397     /* get the pointer to the kernel flavor we need to use */
 398     nb_kernel = select_nbnxn_kernel(nbp->eeltype,
 399                                     nbp->vdwtype,
 400                                     bCalcEner,
 401                                     (plist->haveFreshList && !nb->timers->didPrune[iloc]),
 402                                     nb->dev_info);
 403
 404     /* Kernel launch config:
 405      * - The thread block dimensions match the size of i-clusters, j-clusters,
 406      *   and j-cluster concurrency, in x, y, and z, respectively.
 407      * - The 1D block-grid contains as many blocks as super-clusters.
 408      */
 409     int num_threads_z = 1;
 410     if (nb->dev_info->prop.major == 3 && nb->dev_info->prop.minor == 7)
 411     {
 412         num_threads_z = 2;
 413     }
 414     nblock    = calc_nb_kernel_nblock(plist->nsci, nb->dev_info);
 415
 416     KernelLaunchConfig config;
 417     config.blockSize[0]     = c_clSize;
 418     config.blockSize[1]     = c_clSize;
 419     config.blockSize[2]     = num_threads_z;
 420     config.gridSize[0]      = nblock;
 421     config.sharedMemorySize = calc_shmem_required_nonbonded(num_threads_z, nb->dev_info, nbp);
 422     config.stream           = stream;
 423
 424     if (debug)
 425     {
 426         fprintf(debug, "Non-bonded GPU launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 427                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 428                 "\tShMem: %zu\n",
 429                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 430                 config.gridSize[0], config.gridSize[1], plist->nsci*c_numClPerSupercl,
 431                 c_numClPerSupercl, plist->na_c,
 432                 config.sharedMemorySize);
 433     }
 434
 435     auto      *timingEvent = bDoTime ? t->nb_k[iloc].fetchNextEvent() : nullptr;
 436     const auto kernelArgs  = prepareGpuKernelArguments(nb_kernel, config, adat, nbp, plist, &bCalcFshift);
 437     launchGpuKernel(nb_kernel, config, timingEvent, "k_calc_nb", kernelArgs);
 438
 439     if (bDoTime)
 440     {
 441         t->nb_k[iloc].closeTimingRegion(stream);
 442     }
 443
 444     if (GMX_NATIVE_WINDOWS)
 445     {
 446         /* Windows: force flushing WDDM queue */
 447         cudaStreamQuery(stream);
 448     }
 449 }
 450
 451 /*! Calculates the amount of shared memory required by the CUDA kernel in use. */
 452 static inline int calc_shmem_required_prune(const int num_threads_z)
 453 {
 454     int shmem;
 455
 456     /* i-atom x in shared memory */
 457     shmem  = c_numClPerSupercl * c_clSize * sizeof(float4);
 458     /* cj in shared memory, for each warp separately */
 459     shmem += num_threads_z * c_nbnxnGpuClusterpairSplit * c_nbnxnGpuJgroupSize * sizeof(int);
 460
 461     return shmem;
 462 }
 463
 464 void nbnxn_gpu_launch_kernel_pruneonly(gmx_nbnxn_cuda_t       *nb,
 465                                        int                     iloc,
 466                                        int                     numParts)
 467 {
 468     cu_atomdata_t       *adat    = nb->atdat;
 469     cu_nbparam_t        *nbp     = nb->nbparam;
 470     cu_plist_t          *plist   = nb->plist[iloc];
 471     cu_timers_t         *t       = nb->timers;
 472     cudaStream_t         stream  = nb->stream[iloc];
 473
 474     bool                 bDoTime = nb->bDoTime;
 475
 476     if (plist->haveFreshList)
 477     {
 478         GMX_ASSERT(numParts == 1, "With first pruning we expect 1 part");
 479
 480         /* Set rollingPruningNumParts to signal that it is not set */
 481         plist->rollingPruningNumParts = 0;
 482         plist->rollingPruningPart     = 0;
 483     }
 484     else
 485     {
 486         if (plist->rollingPruningNumParts == 0)
 487         {
 488             plist->rollingPruningNumParts = numParts;
 489         }
 490         else
 491         {
 492             GMX_ASSERT(numParts == plist->rollingPruningNumParts, "It is not allowed to change numParts in between list generation steps");
 493         }
 494     }
 495
 496     /* Use a local variable for part and update in plist, so we can return here
 497      * without duplicating the part increment code.
 498      */
 499     int part = plist->rollingPruningPart;
 500
 501     plist->rollingPruningPart++;
 502     if (plist->rollingPruningPart >= plist->rollingPruningNumParts)
 503     {
 504         plist->rollingPruningPart = 0;
 505     }
 506
 507     /* Compute the number of list entries to prune in this pass */
 508     int numSciInPart = (plist->nsci - part)/numParts;
 509
 510     /* Don't launch the kernel if there is no work to do (not allowed with CUDA) */
 511     if (numSciInPart <= 0)
 512     {
 513         plist->haveFreshList = false;
 514
 515         return;
 516     }
 517
 518     GpuRegionTimer *timer = nullptr;
 519     if (bDoTime)
 520     {
 521         timer = &(plist->haveFreshList ? t->prune_k[iloc] : t->rollingPrune_k[iloc]);
 522     }
 523
 524     /* beginning of timed prune calculation section */
 525     if (bDoTime)
 526     {
 527         timer->openTimingRegion(stream);
 528     }
 529
 530     /* Kernel launch config:
 531      * - The thread block dimensions match the size of i-clusters, j-clusters,
 532      *   and j-cluster concurrency, in x, y, and z, respectively.
 533      * - The 1D block-grid contains as many blocks as super-clusters.
 534      */
 535     int                num_threads_z  = c_cudaPruneKernelJ4Concurrency;
 536     int                nblock         = calc_nb_kernel_nblock(numSciInPart, nb->dev_info);
 537     KernelLaunchConfig config;
 538     config.blockSize[0]     = c_clSize;
 539     config.blockSize[1]     = c_clSize;
 540     config.blockSize[2]     = num_threads_z;
 541     config.gridSize[0]      = nblock;
 542     config.sharedMemorySize = calc_shmem_required_prune(num_threads_z);
 543     config.stream           = stream;
 544
 545     if (debug)
 546     {
 547         fprintf(debug, "Pruning GPU kernel launch configuration:\n\tThread block: %zux%zux%zu\n\t"
 548                 "\tGrid: %zux%zu\n\t#Super-clusters/clusters: %d/%d (%d)\n"
 549                 "\tShMem: %zu\n",
 550                 config.blockSize[0], config.blockSize[1], config.blockSize[2],
 551                 config.gridSize[0], config.gridSize[1], numSciInPart*c_numClPerSupercl,
 552                 c_numClPerSupercl, plist->na_c,
 553                 config.sharedMemorySize);
 554     }
 555
 556     auto          *timingEvent  = bDoTime ? timer->fetchNextEvent() : nullptr;
 557     constexpr char kernelName[] = "k_pruneonly";
 558     const auto    &kernel       = plist->haveFreshList ? nbnxn_kernel_prune_cuda<true> : nbnxn_kernel_prune_cuda<false>;
 559     const auto     kernelArgs   = prepareGpuKernelArguments(kernel, config, adat, nbp, plist, &numParts, &part);
 560     launchGpuKernel(kernel, config, timingEvent, kernelName, kernelArgs);
 561
 562     /* TODO: consider a more elegant way to track which kernel has been called
 563        (combined or separate 1st pass prune, rolling prune). */
 564     if (plist->haveFreshList)
 565     {
 566         plist->haveFreshList         = false;
 567         /* Mark that pruning has been done */
 568         nb->timers->didPrune[iloc] = true;
 569     }
 570     else
 571     {
 572         /* Mark that rolling pruning has been done */
 573         nb->timers->didRollingPrune[iloc] = true;
 574     }
 575
 576     if (bDoTime)
 577     {
 578         timer->closeTimingRegion(stream);
 579     }
 580
 581     if (GMX_NATIVE_WINDOWS)
 582     {
 583         /* Windows: force flushing WDDM queue */
 584         cudaStreamQuery(stream);
 585     }
 586 }
 587
 588 void nbnxn_gpu_launch_cpyback(gmx_nbnxn_cuda_t       *nb,
 589                               const nbnxn_atomdata_t *nbatom,
 590                               int                     flags,
 591                               int                     aloc)
 592 {
 593     cudaError_t stat;
 594     int         adat_begin, adat_len; /* local/nonlocal offset and length used for xq and f */
 595
 596     /* determine interaction locality from atom locality */
 597     int              iloc = gpuAtomToInteractionLocality(aloc);
 598
 599     cu_atomdata_t   *adat    = nb->atdat;
 600     cu_timers_t     *t       = nb->timers;
 601     bool             bDoTime = nb->bDoTime;
 602     cudaStream_t     stream  = nb->stream[iloc];
 603
 604     bool             bCalcEner   = flags & GMX_FORCE_ENERGY;
 605     bool             bCalcFshift = flags & GMX_FORCE_VIRIAL;
 606
 607     /* don't launch non-local copy-back if there was no non-local work to do */
 608     if (canSkipWork(nb, iloc))
 609     {
 610         return;
 611     }
 612
 613     getGpuAtomRange(adat, aloc, adat_begin, adat_len);
 614
 615     /* beginning of timed D2H section */
 616     if (bDoTime)
 617     {
 618         t->nb_d2h[iloc].openTimingRegion(stream);
 619     }
 620
 621     /* With DD the local D2H transfer can only start after the non-local
 622        kernel has finished. */
 623     if (iloc == eintLocal && nb->bUseTwoStreams)
 624     {
 625         stat = cudaStreamWaitEvent(stream, nb->nonlocal_done, 0);
 626         CU_RET_ERR(stat, "cudaStreamWaitEvent on nonlocal_done failed");
 627     }
 628
 629     /* DtoH f */
 630     cu_copy_D2H_async(nbatom->out[0].f + adat_begin * 3, adat->f + adat_begin,
 631                       (adat_len)*sizeof(*adat->f), stream);
 632
 633     /* After the non-local D2H is launched the nonlocal_done event can be
 634        recorded which signals that the local D2H can proceed. This event is not
 635        placed after the non-local kernel because we want the non-local data
 636        back first. */
 637     if (iloc == eintNonlocal)
 638     {
 639         stat = cudaEventRecord(nb->nonlocal_done, stream);
 640         CU_RET_ERR(stat, "cudaEventRecord on nonlocal_done failed");
 641     }
 642
 643     /* only transfer energies in the local stream */
 644     if (LOCAL_I(iloc))
 645     {
 646         /* DtoH fshift */
 647         if (bCalcFshift)
 648         {
 649             cu_copy_D2H_async(nb->nbst.fshift, adat->fshift,
 650                               SHIFTS * sizeof(*nb->nbst.fshift), stream);
 651         }
 652
 653         /* DtoH energies */
 654         if (bCalcEner)
 655         {
 656             cu_copy_D2H_async(nb->nbst.e_lj, adat->e_lj,
 657                               sizeof(*nb->nbst.e_lj), stream);
 658             cu_copy_D2H_async(nb->nbst.e_el, adat->e_el,
 659                               sizeof(*nb->nbst.e_el), stream);
 660         }
 661     }
 662
 663     if (bDoTime)
 664     {
 665         t->nb_d2h[iloc].closeTimingRegion(stream);
 666     }
 667 }
 668
 669 void nbnxn_cuda_set_cacheconfig(const gmx_device_info_t *devinfo)
 670 {
 671     cudaError_t stat;
 672
 673     for (int i = 0; i < eelCuNR; i++)
 674     {
 675         for (int j = 0; j < evdwCuNR; j++)
 676         {
 677             if (devinfo->prop.major >= 3)
 678             {
 679                 /* Default kernel on sm 3.x and later 32/32 kB Shared/L1 */
 680                 cudaFuncSetCacheConfig(nb_kfunc_ener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 681                 cudaFuncSetCacheConfig(nb_kfunc_ener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 682                 cudaFuncSetCacheConfig(nb_kfunc_noener_prune_ptr[i][j], cudaFuncCachePreferEqual);
 683                 stat = cudaFuncSetCacheConfig(nb_kfunc_noener_noprune_ptr[i][j], cudaFuncCachePreferEqual);
 684             }
 685             else
 686             {
 687                 /* On Fermi prefer L1 gives 2% higher performance */
 688                 /* Default kernel on sm_2.x 16/48 kB Shared/L1 */
 689                 cudaFuncSetCacheConfig(nb_kfunc_ener_prune_ptr[i][j], cudaFuncCachePreferL1);
 690                 cudaFuncSetCacheConfig(nb_kfunc_ener_noprune_ptr[i][j], cudaFuncCachePreferL1);
 691                 cudaFuncSetCacheConfig(nb_kfunc_noener_prune_ptr[i][j], cudaFuncCachePreferL1);
 692                 stat = cudaFuncSetCacheConfig(nb_kfunc_noener_noprune_ptr[i][j], cudaFuncCachePreferL1);
 693             }
 694             CU_RET_ERR(stat, "cudaFuncSetCacheConfig failed");
 695         }
 696     }
 697 }