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