src/gromacs/nbnxm/cuda/nbnxm_cuda_data_mgmt.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_data_mgmt.h
  37  *
  38  *  \author Szilard Pall <pall.szilard@gmail.com>
  39  */
  40 #include "gmxpre.h"
  41
  42 #include <assert.h>
  43 #include <stdarg.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46
  47 #include "gromacs/gpu_utils/cudautils.cuh"
  48 #include "gromacs/gpu_utils/gpu_utils.h"
  49 #include "gromacs/gpu_utils/pmalloc_cuda.h"
  50 #include "gromacs/hardware/gpu_hw_info.h"
  51 #include "gromacs/math/vectypes.h"
  52 #include "gromacs/mdlib/force_flags.h"
  53 #include "gromacs/mdtypes/interaction_const.h"
  54 #include "gromacs/mdtypes/md_enums.h"
  55 #include "gromacs/nbnxm/gpu_data_mgmt.h"
  56 #include "gromacs/nbnxm/nbnxm.h"
  57 #include "gromacs/pbcutil/ishift.h"
  58 #include "gromacs/timing/gpu_timing.h"
  59 #include "gromacs/utility/basedefinitions.h"
  60 #include "gromacs/utility/cstringutil.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/real.h"
  63 #include "gromacs/utility/smalloc.h"
  64
  65 #include "nbnxm_cuda.h"
  66 #include "nbnxm_cuda_types.h"
  67
  68 /* This is a heuristically determined parameter for the Kepler
  69  * and Maxwell architectures for the minimum size of ci lists by multiplying
  70  * this constant with the # of multiprocessors on the current device.
  71  * Since the maximum number of blocks per multiprocessor is 16, the ideal
  72  * count for small systems is 32 or 48 blocks per multiprocessor. Because
  73  * there is a bit of fluctuations in the generated block counts, we use
  74  * a target of 44 instead of the ideal value of 48.
  75  */
  76 static unsigned int gpu_min_ci_balanced_factor = 44;
  77
  78 /* Fw. decl. */
  79 static void nbnxn_cuda_clear_e_fshift(gmx_nbnxn_cuda_t *nb);
  80
  81 /* Fw. decl, */
  82 static void nbnxn_cuda_free_nbparam_table(cu_nbparam_t            *nbparam);
  83
  84 /*! \brief Return whether combination rules are used.
  85  *
  86  * \param[in]   pointer to nonbonded paramter struct
  87  * \return      true if combination rules are used in this run, false otherwise
  88  */
  89 static inline bool useLjCombRule(const cu_nbparam_t  *nbparam)
  90 {
  91     return (nbparam->vdwtype == evdwCuCUTCOMBGEOM ||
  92             nbparam->vdwtype == evdwCuCUTCOMBLB);
  93 }
  94
  95 /*! \brief Initialized the Ewald Coulomb correction GPU table.
  96
  97     Tabulates the Ewald Coulomb force and initializes the size/scale
  98     and the table GPU array. If called with an already allocated table,
  99     it just re-uploads the table.
 100  */
 101 static void init_ewald_coulomb_force_table(const interaction_const_t *ic,
 102                                            cu_nbparam_t              *nbp)
 103 {
 104     if (nbp->coulomb_tab != nullptr)
 105     {
 106         nbnxn_cuda_free_nbparam_table(nbp);
 107     }
 108
 109     nbp->coulomb_tab_scale = ic->tabq_scale;
 110     initParamLookupTable(nbp->coulomb_tab, nbp->coulomb_tab_texobj,
 111                          ic->tabq_coul_F, ic->tabq_size);
 112 }
 113
 114
 115 /*! Initializes the atomdata structure first time, it only gets filled at
 116     pair-search. */
 117 static void init_atomdata_first(cu_atomdata_t *ad, int ntypes)
 118 {
 119     cudaError_t stat;
 120
 121     ad->ntypes  = ntypes;
 122     stat        = cudaMalloc((void**)&ad->shift_vec, SHIFTS*sizeof(*ad->shift_vec));
 123     CU_RET_ERR(stat, "cudaMalloc failed on ad->shift_vec");
 124     ad->bShiftVecUploaded = false;
 125
 126     stat = cudaMalloc((void**)&ad->fshift, SHIFTS*sizeof(*ad->fshift));
 127     CU_RET_ERR(stat, "cudaMalloc failed on ad->fshift");
 128
 129     stat = cudaMalloc((void**)&ad->e_lj, sizeof(*ad->e_lj));
 130     CU_RET_ERR(stat, "cudaMalloc failed on ad->e_lj");
 131     stat = cudaMalloc((void**)&ad->e_el, sizeof(*ad->e_el));
 132     CU_RET_ERR(stat, "cudaMalloc failed on ad->e_el");
 133
 134     /* initialize to nullptr poiters to data that is not allocated here and will
 135        need reallocation in nbnxn_cuda_init_atomdata */
 136     ad->xq = nullptr;
 137     ad->f  = nullptr;
 138
 139     /* size -1 indicates that the respective array hasn't been initialized yet */
 140     ad->natoms = -1;
 141     ad->nalloc = -1;
 142 }
 143
 144 /*! Selects the Ewald kernel type, analytical on SM 3.0 and later, tabulated on
 145     earlier GPUs, single or twin cut-off. */
 146 static int pick_ewald_kernel_type(bool                     bTwinCut)
 147 {
 148     bool bUseAnalyticalEwald, bForceAnalyticalEwald, bForceTabulatedEwald;
 149     int  kernel_type;
 150
 151     /* Benchmarking/development environment variables to force the use of
 152        analytical or tabulated Ewald kernel. */
 153     bForceAnalyticalEwald = (getenv("GMX_CUDA_NB_ANA_EWALD") != nullptr);
 154     bForceTabulatedEwald  = (getenv("GMX_CUDA_NB_TAB_EWALD") != nullptr);
 155
 156     if (bForceAnalyticalEwald && bForceTabulatedEwald)
 157     {
 158         gmx_incons("Both analytical and tabulated Ewald CUDA non-bonded kernels "
 159                    "requested through environment variables.");
 160     }
 161
 162     /* By default use analytical Ewald. */
 163     bUseAnalyticalEwald = true;
 164     if (bForceAnalyticalEwald)
 165     {
 166         if (debug)
 167         {
 168             fprintf(debug, "Using analytical Ewald CUDA kernels\n");
 169         }
 170     }
 171     else if (bForceTabulatedEwald)
 172     {
 173         bUseAnalyticalEwald = false;
 174
 175         if (debug)
 176         {
 177             fprintf(debug, "Using tabulated Ewald CUDA kernels\n");
 178         }
 179     }
 180
 181     /* Use twin cut-off kernels if requested by bTwinCut or the env. var.
 182        forces it (use it for debugging/benchmarking only). */
 183     if (!bTwinCut && (getenv("GMX_CUDA_NB_EWALD_TWINCUT") == nullptr))
 184     {
 185         kernel_type = bUseAnalyticalEwald ? eelCuEWALD_ANA : eelCuEWALD_TAB;
 186     }
 187     else
 188     {
 189         kernel_type = bUseAnalyticalEwald ? eelCuEWALD_ANA_TWIN : eelCuEWALD_TAB_TWIN;
 190     }
 191
 192     return kernel_type;
 193 }
 194
 195 /*! Copies all parameters related to the cut-off from ic to nbp */
 196 static void set_cutoff_parameters(cu_nbparam_t              *nbp,
 197                                   const interaction_const_t *ic,
 198                                   const NbnxnListParameters *listParams)
 199 {
 200     nbp->ewald_beta        = ic->ewaldcoeff_q;
 201     nbp->sh_ewald          = ic->sh_ewald;
 202     nbp->epsfac            = ic->epsfac;
 203     nbp->two_k_rf          = 2.0 * ic->k_rf;
 204     nbp->c_rf              = ic->c_rf;
 205     nbp->rvdw_sq           = ic->rvdw * ic->rvdw;
 206     nbp->rcoulomb_sq       = ic->rcoulomb * ic->rcoulomb;
 207     nbp->rlistOuter_sq     = listParams->rlistOuter * listParams->rlistOuter;
 208     nbp->rlistInner_sq     = listParams->rlistInner * listParams->rlistInner;
 209     nbp->useDynamicPruning = listParams->useDynamicPruning;
 210
 211     nbp->sh_lj_ewald       = ic->sh_lj_ewald;
 212     nbp->ewaldcoeff_lj     = ic->ewaldcoeff_lj;
 213
 214     nbp->rvdw_switch       = ic->rvdw_switch;
 215     nbp->dispersion_shift  = ic->dispersion_shift;
 216     nbp->repulsion_shift   = ic->repulsion_shift;
 217     nbp->vdw_switch        = ic->vdw_switch;
 218 }
 219
 220 /*! Initializes the nonbonded parameter data structure. */
 221 static void init_nbparam(cu_nbparam_t                   *nbp,
 222                          const interaction_const_t      *ic,
 223                          const NbnxnListParameters      *listParams,
 224                          const nbnxn_atomdata_t::Params &nbatParams)
 225 {
 226     int         ntypes;
 227
 228     ntypes  = nbatParams.numTypes;
 229
 230     set_cutoff_parameters(nbp, ic, listParams);
 231
 232     /* The kernel code supports LJ combination rules (geometric and LB) for
 233      * all kernel types, but we only generate useful combination rule kernels.
 234      * We currently only use LJ combination rule (geometric and LB) kernels
 235      * for plain cut-off LJ. On Maxwell the force only kernels speed up 15%
 236      * with PME and 20% with RF, the other kernels speed up about half as much.
 237      * For LJ force-switch the geometric rule would give 7% speed-up, but this
 238      * combination is rarely used. LJ force-switch with LB rule is more common,
 239      * but gives only 1% speed-up.
 240      */
 241     if (ic->vdwtype == evdwCUT)
 242     {
 243         switch (ic->vdw_modifier)
 244         {
 245             case eintmodNONE:
 246             case eintmodPOTSHIFT:
 247                 switch (nbatParams.comb_rule)
 248                 {
 249                     case ljcrNONE:
 250                         nbp->vdwtype = evdwCuCUT;
 251                         break;
 252                     case ljcrGEOM:
 253                         nbp->vdwtype = evdwCuCUTCOMBGEOM;
 254                         break;
 255                     case ljcrLB:
 256                         nbp->vdwtype = evdwCuCUTCOMBLB;
 257                         break;
 258                     default:
 259                         gmx_incons("The requested LJ combination rule is not implemented in the CUDA GPU accelerated kernels!");
 260                 }
 261                 break;
 262             case eintmodFORCESWITCH:
 263                 nbp->vdwtype = evdwCuFSWITCH;
 264                 break;
 265             case eintmodPOTSWITCH:
 266                 nbp->vdwtype = evdwCuPSWITCH;
 267                 break;
 268             default:
 269                 gmx_incons("The requested VdW interaction modifier is not implemented in the CUDA GPU accelerated kernels!");
 270         }
 271     }
 272     else if (ic->vdwtype == evdwPME)
 273     {
 274         if (ic->ljpme_comb_rule == ljcrGEOM)
 275         {
 276             assert(nbatParams.comb_rule == ljcrGEOM);
 277             nbp->vdwtype = evdwCuEWALDGEOM;
 278         }
 279         else
 280         {
 281             assert(nbatParams.comb_rule == ljcrLB);
 282             nbp->vdwtype = evdwCuEWALDLB;
 283         }
 284     }
 285     else
 286     {
 287         gmx_incons("The requested VdW type is not implemented in the CUDA GPU accelerated kernels!");
 288     }
 289
 290     if (ic->eeltype == eelCUT)
 291     {
 292         nbp->eeltype = eelCuCUT;
 293     }
 294     else if (EEL_RF(ic->eeltype))
 295     {
 296         nbp->eeltype = eelCuRF;
 297     }
 298     else if ((EEL_PME(ic->eeltype) || ic->eeltype == eelEWALD))
 299     {
 300         /* Initially rcoulomb == rvdw, so it's surely not twin cut-off. */
 301         nbp->eeltype = pick_ewald_kernel_type(false);
 302     }
 303     else
 304     {
 305         /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
 306         gmx_incons("The requested electrostatics type is not implemented in the CUDA GPU accelerated kernels!");
 307     }
 308
 309     /* generate table for PME */
 310     nbp->coulomb_tab = nullptr;
 311     if (nbp->eeltype == eelCuEWALD_TAB || nbp->eeltype == eelCuEWALD_TAB_TWIN)
 312     {
 313         init_ewald_coulomb_force_table(ic, nbp);
 314     }
 315
 316     /* set up LJ parameter lookup table */
 317     if (!useLjCombRule(nbp))
 318     {
 319         initParamLookupTable(nbp->nbfp, nbp->nbfp_texobj,
 320                              nbatParams.nbfp.data(), 2*ntypes*ntypes);
 321     }
 322
 323     /* set up LJ-PME parameter lookup table */
 324     if (ic->vdwtype == evdwPME)
 325     {
 326         initParamLookupTable(nbp->nbfp_comb, nbp->nbfp_comb_texobj,
 327                              nbatParams.nbfp_comb.data(), 2*ntypes);
 328     }
 329 }
 330
 331 /*! Re-generate the GPU Ewald force table, resets rlist, and update the
 332  *  electrostatic type switching to twin cut-off (or back) if needed. */
 333 void nbnxn_gpu_pme_loadbal_update_param(const nonbonded_verlet_t    *nbv,
 334                                         const interaction_const_t   *ic,
 335                                         const NbnxnListParameters   *listParams)
 336 {
 337     if (!nbv || nbv->grp[0].kernel_type != nbnxnk8x8x8_GPU)
 338     {
 339         return;
 340     }
 341     gmx_nbnxn_cuda_t *nb    = nbv->gpu_nbv;
 342     cu_nbparam_t     *nbp   = nb->nbparam;
 343
 344     set_cutoff_parameters(nbp, ic, listParams);
 345
 346     nbp->eeltype        = pick_ewald_kernel_type(ic->rcoulomb != ic->rvdw);
 347
 348     init_ewald_coulomb_force_table(ic, nb->nbparam);
 349 }
 350
 351 /*! Initializes the pair list data structure. */
 352 static void init_plist(cu_plist_t *pl)
 353 {
 354     /* initialize to nullptr pointers to data that is not allocated here and will
 355        need reallocation in nbnxn_gpu_init_pairlist */
 356     pl->sci      = nullptr;
 357     pl->cj4      = nullptr;
 358     pl->imask    = nullptr;
 359     pl->excl     = nullptr;
 360
 361     /* size -1 indicates that the respective array hasn't been initialized yet */
 362     pl->na_c           = -1;
 363     pl->nsci           = -1;
 364     pl->sci_nalloc     = -1;
 365     pl->ncj4           = -1;
 366     pl->cj4_nalloc     = -1;
 367     pl->nimask         = -1;
 368     pl->imask_nalloc   = -1;
 369     pl->nexcl          = -1;
 370     pl->excl_nalloc    = -1;
 371     pl->haveFreshList  = false;
 372 }
 373
 374 /*! Initializes the timer data structure. */
 375 static void init_timers(cu_timers_t *t, bool bUseTwoStreams)
 376 {
 377     /* The non-local counters/stream (second in the array) are needed only with DD. */
 378     for (int i = 0; i <= (bUseTwoStreams ? 1 : 0); i++)
 379     {
 380         t->didPairlistH2D[i]  = false;
 381         t->didPrune[i]        = false;
 382         t->didRollingPrune[i] = false;
 383     }
 384 }
 385
 386 /*! Initializes the timings data structure. */
 387 static void init_timings(gmx_wallclock_gpu_nbnxn_t *t)
 388 {
 389     int i, j;
 390
 391     t->nb_h2d_t = 0.0;
 392     t->nb_d2h_t = 0.0;
 393     t->nb_c     = 0;
 394     t->pl_h2d_t = 0.0;
 395     t->pl_h2d_c = 0;
 396     for (i = 0; i < 2; i++)
 397     {
 398         for (j = 0; j < 2; j++)
 399         {
 400             t->ktime[i][j].t = 0.0;
 401             t->ktime[i][j].c = 0;
 402         }
 403     }
 404     t->pruneTime.c        = 0;
 405     t->pruneTime.t        = 0.0;
 406     t->dynamicPruneTime.c = 0;
 407     t->dynamicPruneTime.t = 0.0;
 408 }
 409
 410 /*! Initializes simulation constant data. */
 411 static void nbnxn_cuda_init_const(gmx_nbnxn_cuda_t               *nb,
 412                                   const interaction_const_t      *ic,
 413                                   const NbnxnListParameters      *listParams,
 414                                   const nbnxn_atomdata_t::Params &nbatParams)
 415 {
 416     init_atomdata_first(nb->atdat, nbatParams.numTypes);
 417     init_nbparam(nb->nbparam, ic, listParams, nbatParams);
 418
 419     /* clear energy and shift force outputs */
 420     nbnxn_cuda_clear_e_fshift(nb);
 421 }
 422
 423 void nbnxn_gpu_init(gmx_nbnxn_cuda_t         **p_nb,
 424                     const gmx_device_info_t   *deviceInfo,
 425                     const interaction_const_t *ic,
 426                     const NbnxnListParameters *listParams,
 427                     const nbnxn_atomdata_t    *nbat,
 428                     int                        /*rank*/,
 429                     gmx_bool                   bLocalAndNonlocal)
 430 {
 431     cudaError_t       stat;
 432     gmx_nbnxn_cuda_t *nb;
 433
 434     if (p_nb == nullptr)
 435     {
 436         return;
 437     }
 438
 439     snew(nb, 1);
 440     snew(nb->atdat, 1);
 441     snew(nb->nbparam, 1);
 442     snew(nb->plist[eintLocal], 1);
 443     if (bLocalAndNonlocal)
 444     {
 445         snew(nb->plist[eintNonlocal], 1);
 446     }
 447
 448     nb->bUseTwoStreams = bLocalAndNonlocal;
 449
 450     nb->timers = new cu_timers_t();
 451     snew(nb->timings, 1);
 452
 453     /* init nbst */
 454     pmalloc((void**)&nb->nbst.e_lj, sizeof(*nb->nbst.e_lj));
 455     pmalloc((void**)&nb->nbst.e_el, sizeof(*nb->nbst.e_el));
 456     pmalloc((void**)&nb->nbst.fshift, SHIFTS * sizeof(*nb->nbst.fshift));
 457
 458     init_plist(nb->plist[eintLocal]);
 459
 460     /* set device info, just point it to the right GPU among the detected ones */
 461     nb->dev_info = deviceInfo;
 462
 463     /* local/non-local GPU streams */
 464     stat = cudaStreamCreate(&nb->stream[eintLocal]);
 465     CU_RET_ERR(stat, "cudaStreamCreate on stream[eintLocal] failed");
 466     if (nb->bUseTwoStreams)
 467     {
 468         init_plist(nb->plist[eintNonlocal]);
 469
 470         /* Note that the device we're running on does not have to support
 471          * priorities, because we are querying the priority range which in this
 472          * case will be a single value.
 473          */
 474         int highest_priority;
 475         stat = cudaDeviceGetStreamPriorityRange(nullptr, &highest_priority);
 476         CU_RET_ERR(stat, "cudaDeviceGetStreamPriorityRange failed");
 477
 478         stat = cudaStreamCreateWithPriority(&nb->stream[eintNonlocal],
 479                                             cudaStreamDefault,
 480                                             highest_priority);
 481         CU_RET_ERR(stat, "cudaStreamCreateWithPriority on stream[eintNonlocal] failed");
 482     }
 483
 484     /* init events for sychronization (timing disabled for performance reasons!) */
 485     stat = cudaEventCreateWithFlags(&nb->nonlocal_done, cudaEventDisableTiming);
 486     CU_RET_ERR(stat, "cudaEventCreate on nonlocal_done failed");
 487     stat = cudaEventCreateWithFlags(&nb->misc_ops_and_local_H2D_done, cudaEventDisableTiming);
 488     CU_RET_ERR(stat, "cudaEventCreate on misc_ops_and_local_H2D_done failed");
 489
 490     /* WARNING: CUDA timings are incorrect with multiple streams.
 491      *          This is the main reason why they are disabled by default.
 492      */
 493     // TODO: Consider turning on by default when we can detect nr of streams.
 494     nb->bDoTime = (getenv("GMX_ENABLE_GPU_TIMING") != nullptr);
 495
 496     if (nb->bDoTime)
 497     {
 498         init_timers(nb->timers, nb->bUseTwoStreams);
 499         init_timings(nb->timings);
 500     }
 501
 502     /* set the kernel type for the current GPU */
 503     /* pick L1 cache configuration */
 504     nbnxn_cuda_set_cacheconfig();
 505
 506     nbnxn_cuda_init_const(nb, ic, listParams, nbat->params());
 507
 508     *p_nb = nb;
 509
 510     if (debug)
 511     {
 512         fprintf(debug, "Initialized CUDA data structures.\n");
 513     }
 514 }
 515
 516 void nbnxn_gpu_init_pairlist(gmx_nbnxn_cuda_t       *nb,
 517                              const NbnxnPairlistGpu *h_plist,
 518                              int                     iloc)
 519 {
 520     char          sbuf[STRLEN];
 521     bool          bDoTime    =  (nb->bDoTime && !h_plist->sci.empty());
 522     cudaStream_t  stream     = nb->stream[iloc];
 523     cu_plist_t   *d_plist    = nb->plist[iloc];
 524
 525     if (d_plist->na_c < 0)
 526     {
 527         d_plist->na_c = h_plist->na_ci;
 528     }
 529     else
 530     {
 531         if (d_plist->na_c != h_plist->na_ci)
 532         {
 533             sprintf(sbuf, "In cu_init_plist: the #atoms per cell has changed (from %d to %d)",
 534                     d_plist->na_c, h_plist->na_ci);
 535             gmx_incons(sbuf);
 536         }
 537     }
 538
 539     if (bDoTime)
 540     {
 541         nb->timers->pl_h2d[iloc].openTimingRegion(stream);
 542         nb->timers->didPairlistH2D[iloc] = true;
 543     }
 544
 545     Context context = nullptr;
 546
 547     reallocateDeviceBuffer(&d_plist->sci, h_plist->sci.size(),
 548                            &d_plist->nsci, &d_plist->sci_nalloc, context);
 549     copyToDeviceBuffer(&d_plist->sci, h_plist->sci.data(), 0, h_plist->sci.size(),
 550                        stream, GpuApiCallBehavior::Async,
 551                        bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
 552
 553     reallocateDeviceBuffer(&d_plist->cj4, h_plist->cj4.size(),
 554                            &d_plist->ncj4, &d_plist->cj4_nalloc, context);
 555     copyToDeviceBuffer(&d_plist->cj4, h_plist->cj4.data(), 0, h_plist->cj4.size(),
 556                        stream, GpuApiCallBehavior::Async,
 557                        bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
 558
 559     reallocateDeviceBuffer(&d_plist->imask, h_plist->cj4.size()*c_nbnxnGpuClusterpairSplit,
 560                            &d_plist->nimask, &d_plist->imask_nalloc, context);
 561
 562     reallocateDeviceBuffer(&d_plist->excl, h_plist->excl.size(),
 563                            &d_plist->nexcl, &d_plist->excl_nalloc, context);
 564     copyToDeviceBuffer(&d_plist->excl, h_plist->excl.data(), 0, h_plist->excl.size(),
 565                        stream, GpuApiCallBehavior::Async,
 566                        bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
 567
 568     if (bDoTime)
 569     {
 570         nb->timers->pl_h2d[iloc].closeTimingRegion(stream);
 571     }
 572
 573     /* the next use of thist list we be the first one, so we need to prune */
 574     d_plist->haveFreshList = true;
 575 }
 576
 577 void nbnxn_gpu_upload_shiftvec(gmx_nbnxn_cuda_t       *nb,
 578                                const nbnxn_atomdata_t *nbatom)
 579 {
 580     cu_atomdata_t *adat  = nb->atdat;
 581     cudaStream_t   ls    = nb->stream[eintLocal];
 582
 583     /* only if we have a dynamic box */
 584     if (nbatom->bDynamicBox || !adat->bShiftVecUploaded)
 585     {
 586         cu_copy_H2D_async(adat->shift_vec, nbatom->shift_vec.data(),
 587                           SHIFTS * sizeof(*adat->shift_vec), ls);
 588         adat->bShiftVecUploaded = true;
 589     }
 590 }
 591
 592 /*! Clears the first natoms_clear elements of the GPU nonbonded force output array. */
 593 static void nbnxn_cuda_clear_f(gmx_nbnxn_cuda_t *nb, int natoms_clear)
 594 {
 595     cudaError_t    stat;
 596     cu_atomdata_t *adat  = nb->atdat;
 597     cudaStream_t   ls    = nb->stream[eintLocal];
 598
 599     stat = cudaMemsetAsync(adat->f, 0, natoms_clear * sizeof(*adat->f), ls);
 600     CU_RET_ERR(stat, "cudaMemsetAsync on f falied");
 601 }
 602
 603 /*! Clears nonbonded shift force output array and energy outputs on the GPU. */
 604 static void nbnxn_cuda_clear_e_fshift(gmx_nbnxn_cuda_t *nb)
 605 {
 606     cudaError_t    stat;
 607     cu_atomdata_t *adat  = nb->atdat;
 608     cudaStream_t   ls    = nb->stream[eintLocal];
 609
 610     stat = cudaMemsetAsync(adat->fshift, 0, SHIFTS * sizeof(*adat->fshift), ls);
 611     CU_RET_ERR(stat, "cudaMemsetAsync on fshift falied");
 612     stat = cudaMemsetAsync(adat->e_lj, 0, sizeof(*adat->e_lj), ls);
 613     CU_RET_ERR(stat, "cudaMemsetAsync on e_lj falied");
 614     stat = cudaMemsetAsync(adat->e_el, 0, sizeof(*adat->e_el), ls);
 615     CU_RET_ERR(stat, "cudaMemsetAsync on e_el falied");
 616 }
 617
 618 void nbnxn_gpu_clear_outputs(gmx_nbnxn_cuda_t *nb, int flags)
 619 {
 620     nbnxn_cuda_clear_f(nb, nb->atdat->natoms);
 621     /* clear shift force array and energies if the outputs were
 622        used in the current step */
 623     if (flags & GMX_FORCE_VIRIAL)
 624     {
 625         nbnxn_cuda_clear_e_fshift(nb);
 626     }
 627 }
 628
 629 void nbnxn_gpu_init_atomdata(gmx_nbnxn_cuda_t              *nb,
 630                              const nbnxn_atomdata_t        *nbat)
 631 {
 632     cudaError_t    stat;
 633     int            nalloc, natoms;
 634     bool           realloced;
 635     bool           bDoTime   = nb->bDoTime;
 636     cu_timers_t   *timers    = nb->timers;
 637     cu_atomdata_t *d_atdat   = nb->atdat;
 638     cudaStream_t   ls        = nb->stream[eintLocal];
 639
 640     natoms    = nbat->numAtoms();
 641     realloced = false;
 642
 643     if (bDoTime)
 644     {
 645         /* time async copy */
 646         timers->atdat.openTimingRegion(ls);
 647     }
 648
 649     /* need to reallocate if we have to copy more atoms than the amount of space
 650        available and only allocate if we haven't initialized yet, i.e d_atdat->natoms == -1 */
 651     if (natoms > d_atdat->nalloc)
 652     {
 653         nalloc = over_alloc_small(natoms);
 654
 655         /* free up first if the arrays have already been initialized */
 656         if (d_atdat->nalloc != -1)
 657         {
 658             freeDeviceBuffer(&d_atdat->f);
 659             freeDeviceBuffer(&d_atdat->xq);
 660             freeDeviceBuffer(&d_atdat->atom_types);
 661             freeDeviceBuffer(&d_atdat->lj_comb);
 662         }
 663
 664         stat = cudaMalloc((void **)&d_atdat->f, nalloc*sizeof(*d_atdat->f));
 665         CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->f");
 666         stat = cudaMalloc((void **)&d_atdat->xq, nalloc*sizeof(*d_atdat->xq));
 667         CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->xq");
 668         if (useLjCombRule(nb->nbparam))
 669         {
 670             stat = cudaMalloc((void **)&d_atdat->lj_comb, nalloc*sizeof(*d_atdat->lj_comb));
 671             CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->lj_comb");
 672         }
 673         else
 674         {
 675             stat = cudaMalloc((void **)&d_atdat->atom_types, nalloc*sizeof(*d_atdat->atom_types));
 676             CU_RET_ERR(stat, "cudaMalloc failed on d_atdat->atom_types");
 677         }
 678
 679         d_atdat->nalloc = nalloc;
 680         realloced       = true;
 681     }
 682
 683     d_atdat->natoms       = natoms;
 684     d_atdat->natoms_local = nbat->natoms_local;
 685
 686     /* need to clear GPU f output if realloc happened */
 687     if (realloced)
 688     {
 689         nbnxn_cuda_clear_f(nb, nalloc);
 690     }
 691
 692     if (useLjCombRule(nb->nbparam))
 693     {
 694         cu_copy_H2D_async(d_atdat->lj_comb, nbat->params().lj_comb.data(),
 695                           natoms*sizeof(*d_atdat->lj_comb), ls);
 696     }
 697     else
 698     {
 699         cu_copy_H2D_async(d_atdat->atom_types, nbat->params().type.data(),
 700                           natoms*sizeof(*d_atdat->atom_types), ls);
 701     }
 702
 703     if (bDoTime)
 704     {
 705         timers->atdat.closeTimingRegion(ls);
 706     }
 707 }
 708
 709 static void nbnxn_cuda_free_nbparam_table(cu_nbparam_t            *nbparam)
 710 {
 711     if (nbparam->eeltype == eelCuEWALD_TAB || nbparam->eeltype == eelCuEWALD_TAB_TWIN)
 712     {
 713         destroyParamLookupTable(nbparam->coulomb_tab, nbparam->coulomb_tab_texobj);
 714     }
 715 }
 716
 717 void nbnxn_gpu_free(gmx_nbnxn_cuda_t *nb)
 718 {
 719     cudaError_t      stat;
 720     cu_atomdata_t   *atdat;
 721     cu_nbparam_t    *nbparam;
 722
 723     if (nb == nullptr)
 724     {
 725         return;
 726     }
 727
 728     atdat       = nb->atdat;
 729     nbparam     = nb->nbparam;
 730
 731     nbnxn_cuda_free_nbparam_table(nbparam);
 732
 733     stat = cudaEventDestroy(nb->nonlocal_done);
 734     CU_RET_ERR(stat, "cudaEventDestroy failed on timers->nonlocal_done");
 735     stat = cudaEventDestroy(nb->misc_ops_and_local_H2D_done);
 736     CU_RET_ERR(stat, "cudaEventDestroy failed on timers->misc_ops_and_local_H2D_done");
 737
 738     delete nb->timers;
 739     if (nb->bDoTime)
 740     {
 741         /* The non-local counters/stream (second in the array) are needed only with DD. */
 742         for (int i = 0; i <= (nb->bUseTwoStreams ? 1 : 0); i++)
 743         {
 744             stat = cudaStreamDestroy(nb->stream[i]);
 745             CU_RET_ERR(stat, "cudaStreamDestroy failed on stream");
 746         }
 747     }
 748
 749     if (!useLjCombRule(nb->nbparam))
 750     {
 751         destroyParamLookupTable(nbparam->nbfp, nbparam->nbfp_texobj);
 752
 753     }
 754
 755     if (nbparam->vdwtype == evdwCuEWALDGEOM || nbparam->vdwtype == evdwCuEWALDLB)
 756     {
 757         destroyParamLookupTable(nbparam->nbfp_comb, nbparam->nbfp_comb_texobj);
 758     }
 759
 760     stat = cudaFree(atdat->shift_vec);
 761     CU_RET_ERR(stat, "cudaFree failed on atdat->shift_vec");
 762     stat = cudaFree(atdat->fshift);
 763     CU_RET_ERR(stat, "cudaFree failed on atdat->fshift");
 764
 765     stat = cudaFree(atdat->e_lj);
 766     CU_RET_ERR(stat, "cudaFree failed on atdat->e_lj");
 767     stat = cudaFree(atdat->e_el);
 768     CU_RET_ERR(stat, "cudaFree failed on atdat->e_el");
 769
 770     freeDeviceBuffer(&atdat->f);
 771     freeDeviceBuffer(&atdat->xq);
 772     freeDeviceBuffer(&atdat->atom_types);
 773     freeDeviceBuffer(&atdat->lj_comb);
 774
 775     /* Free plist */
 776     auto *plist = nb->plist[eintLocal];
 777     freeDeviceBuffer(&plist->sci);
 778     freeDeviceBuffer(&plist->cj4);
 779     freeDeviceBuffer(&plist->imask);
 780     freeDeviceBuffer(&plist->excl);
 781     sfree(plist);
 782     if (nb->bUseTwoStreams)
 783     {
 784         auto *plist_nl = nb->plist[eintNonlocal];
 785         freeDeviceBuffer(&plist_nl->sci);
 786         freeDeviceBuffer(&plist_nl->cj4);
 787         freeDeviceBuffer(&plist_nl->imask);
 788         freeDeviceBuffer(&plist_nl->excl);
 789         sfree(plist_nl);
 790     }
 791
 792     /* Free nbst */
 793     pfree(nb->nbst.e_lj);
 794     nb->nbst.e_lj = nullptr;
 795
 796     pfree(nb->nbst.e_el);
 797     nb->nbst.e_el = nullptr;
 798
 799     pfree(nb->nbst.fshift);
 800     nb->nbst.fshift = nullptr;
 801
 802     sfree(atdat);
 803     sfree(nbparam);
 804     sfree(nb->timings);
 805     sfree(nb);
 806
 807     if (debug)
 808     {
 809         fprintf(debug, "Cleaned up CUDA data structures.\n");
 810     }
 811 }
 812
 813 //! This function is documented in the header file
 814 gmx_wallclock_gpu_nbnxn_t *nbnxn_gpu_get_timings(gmx_nbnxn_cuda_t *nb)
 815 {
 816     return (nb != nullptr && nb->bDoTime) ? nb->timings : nullptr;
 817 }
 818
 819 void nbnxn_gpu_reset_timings(nonbonded_verlet_t* nbv)
 820 {
 821     if (nbv->gpu_nbv && nbv->gpu_nbv->bDoTime)
 822     {
 823         init_timings(nbv->gpu_nbv->timings);
 824     }
 825 }
 826
 827 int nbnxn_gpu_min_ci_balanced(gmx_nbnxn_cuda_t *nb)
 828 {
 829     return nb != nullptr ?
 830            gpu_min_ci_balanced_factor*nb->dev_info->prop.multiProcessorCount : 0;
 831
 832 }
 833
 834 gmx_bool nbnxn_gpu_is_kernel_ewald_analytical(const gmx_nbnxn_cuda_t *nb)
 835 {
 836     return ((nb->nbparam->eeltype == eelCuEWALD_ANA) ||
 837             (nb->nbparam->eeltype == eelCuEWALD_ANA_TWIN));
 838 }
 839
 840 void *nbnxn_gpu_get_command_stream(gmx_nbnxn_gpu_t *nb,
 841                                    int              iloc)
 842 {
 843     assert(nb);
 844
 845     return static_cast<void *>(&nb->stream[iloc]);
 846 }
 847
 848 void *nbnxn_gpu_get_xq(gmx_nbnxn_gpu_t *nb)
 849 {
 850     assert(nb);
 851
 852     return static_cast<void *>(nb->atdat->xq);
 853 }
 854
 855 void *nbnxn_gpu_get_f(gmx_nbnxn_gpu_t *nb)
 856 {
 857     assert(nb);
 858
 859     return static_cast<void *>(nb->atdat->f);
 860 }
 861
 862 rvec *nbnxn_gpu_get_fshift(gmx_nbnxn_gpu_t *nb)
 863 {
 864     assert(nb);
 865
 866     return reinterpret_cast<rvec *>(nb->atdat->fshift);
 867 }