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