src/gromacs/mdlib/nbnxn_ocl/nbnxn_ocl_data_mgmt.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  *  \brief Define OpenCL implementation of nbnxn_gpu_data_mgmt.h
  37  *
  38  *  \author Anca Hamuraru <anca@streamcomputing.eu>
  39  *  \author Dimitrios Karkoulis <dimitris.karkoulis@gmail.com>
  40  *  \author Teemu Virolainen <teemu@streamcomputing.eu>
  41  *  \author Szilárd Páll <pall.szilard@gmail.com>
  42  */
  43 #include "gmxpre.h"
  44
  45 #include <assert.h>
  46 #include <stdarg.h>
  47 #include <stdio.h>
  48 #include <stdlib.h>
  49 #include <string.h>
  50
  51 #include <cmath>
  52
  53 #include "gromacs/gpu_utils/gpu_utils.h"
  54 #include "gromacs/gpu_utils/oclutils.h"
  55 #include "gromacs/hardware/gpu_hw_info.h"
  56 #include "gromacs/math/vectypes.h"
  57 #include "gromacs/mdlib/force_flags.h"
  58 #include "gromacs/mdlib/nb_verlet.h"
  59 #include "gromacs/mdlib/nbnxn_consts.h"
  60 #include "gromacs/mdlib/nbnxn_gpu.h"
  61 #include "gromacs/mdlib/nbnxn_gpu_common_utils.h"
  62 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
  63 #include "gromacs/mdlib/nbnxn_gpu_jit_support.h"
  64 #include "gromacs/mdtypes/interaction_const.h"
  65 #include "gromacs/mdtypes/md_enums.h"
  66 #include "gromacs/pbcutil/ishift.h"
  67 #include "gromacs/timing/gpu_timing.h"
  68 #include "gromacs/utility/cstringutil.h"
  69 #include "gromacs/utility/fatalerror.h"
  70 #include "gromacs/utility/gmxassert.h"
  71 #include "gromacs/utility/real.h"
  72 #include "gromacs/utility/smalloc.h"
  73
  74 #include "nbnxn_ocl_internal.h"
  75 #include "nbnxn_ocl_types.h"
  76
  77 /*! \brief This parameter should be determined heuristically from the
  78  * kernel execution times
  79  *
  80  * This value is best for small systems on a single AMD Radeon R9 290X
  81  * (and about 5% faster than 40, which is the default for CUDA
  82  * devices). Larger simulation systems were quite insensitive to the
  83  * value of this parameter.
  84  */
  85 static unsigned int gpu_min_ci_balanced_factor = 50;
  86
  87
  88 /*! \brief Returns true if LJ combination rules are used in the non-bonded kernels.
  89  *
  90  * Full doc in nbnxn_ocl_internal.h */
  91 bool useLjCombRule(int vdwType)
  92 {
  93     return (vdwType == evdwOclCUTCOMBGEOM ||
  94             vdwType == evdwOclCUTCOMBLB);
  95 }
  96
  97 /*! \brief Reallocation device buffers
  98  *
  99  *  Reallocation of the memory pointed by d_ptr and copying of the data from
 100  *  the location pointed by h_src host-side pointer is done. Allocation is
 101  *  buffered and therefore freeing is only needed if the previously allocated
 102  *  space is not enough.
 103  *  The H2D copy is launched in command queue s and can be done synchronously or
 104  *  asynchronously (the default is the latter).
 105  *  If copy_event is not NULL, on return it will contain an event object
 106  *  identifying the H2D copy. The event can further be used to queue a wait
 107  *  for this operation or to query profiling information.
 108  *  OpenCL equivalent of cu_realloc_buffered.
 109  */
 110 static void ocl_realloc_buffered(cl_mem *d_dest, void *h_src,
 111                                  size_t type_size,
 112                                  int *curr_size, int *curr_alloc_size,
 113                                  int req_size,
 114                                  cl_context context,
 115                                  cl_command_queue s,
 116                                  bool bAsync = true,
 117                                  cl_event *copy_event = NULL)
 118 {
 119     if (d_dest == NULL || req_size < 0)
 120     {
 121         return;
 122     }
 123
 124     /* reallocate only if the data does not fit = allocation size is smaller
 125        than the current requested size */
 126     if (req_size > *curr_alloc_size)
 127     {
 128         cl_int gmx_unused cl_error;
 129
 130         /* only free if the array has already been initialized */
 131         if (*curr_alloc_size >= 0)
 132         {
 133             freeDeviceBuffer(d_dest);
 134         }
 135
 136         *curr_alloc_size = over_alloc_large(req_size);
 137
 138         *d_dest = clCreateBuffer(context, CL_MEM_READ_WRITE, *curr_alloc_size * type_size, NULL, &cl_error);
 139         assert(cl_error == CL_SUCCESS);
 140         // TODO: handle errors, check clCreateBuffer flags
 141     }
 142
 143     /* size could have changed without actual reallocation */
 144     *curr_size = req_size;
 145
 146     /* upload to device */
 147     if (h_src)
 148     {
 149         if (bAsync)
 150         {
 151             ocl_copy_H2D_async(*d_dest, h_src, 0, *curr_size * type_size, s, copy_event);
 152         }
 153         else
 154         {
 155             ocl_copy_H2D_sync(*d_dest, h_src,  0, *curr_size * type_size, s);
 156         }
 157     }
 158 }
 159
 160 /*! \brief Tabulates the Ewald Coulomb force and initializes the size/scale
 161  * and the table GPU array.
 162  *
 163  * If called with an already allocated table, it just re-uploads the
 164  * table.
 165  */
 166 static void init_ewald_coulomb_force_table(const interaction_const_t       *ic,
 167                                            cl_nbparam_t                    *nbp,
 168                                            const gmx_device_runtime_data_t *runData)
 169 {
 170     cl_mem       coul_tab;
 171
 172     cl_int       cl_error;
 173
 174     if (nbp->coulomb_tab_climg2d != NULL)
 175     {
 176         freeDeviceBuffer(&(nbp->coulomb_tab_climg2d));
 177     }
 178
 179     /* Switched from using textures to using buffers */
 180     // TODO: decide which alternative is most efficient - textures or buffers.
 181     /*
 182        cl_image_format array_format;
 183
 184        array_format.image_channel_data_type = CL_FLOAT;
 185        array_format.image_channel_order     = CL_R;
 186
 187        coul_tab = clCreateImage2D(runData->context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
 188        &array_format, tabsize, 1, 0, ftmp, &cl_error);
 189      */
 190
 191     coul_tab = clCreateBuffer(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, ic->tabq_size*sizeof(cl_float), ic->tabq_coul_F, &cl_error);
 192     assert(cl_error == CL_SUCCESS);
 193     // TODO: handle errors, check clCreateBuffer flags
 194
 195     nbp->coulomb_tab_climg2d  = coul_tab;
 196     nbp->coulomb_tab_scale    = ic->tabq_scale;
 197 }
 198
 199
 200 /*! \brief Initializes the atomdata structure first time, it only gets filled at
 201     pair-search.
 202  */
 203 static void init_atomdata_first(cl_atomdata_t *ad, int ntypes, gmx_device_runtime_data_t *runData)
 204 {
 205     cl_int cl_error;
 206
 207     ad->ntypes  = ntypes;
 208
 209     /* An element of the shift_vec device buffer has the same size as one element
 210        of the host side shift_vec buffer. */
 211     ad->shift_vec_elem_size = sizeof(*(((nbnxn_atomdata_t*)0)->shift_vec));
 212
 213     // TODO: handle errors, check clCreateBuffer flags
 214     ad->shift_vec = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, SHIFTS * ad->shift_vec_elem_size, NULL, &cl_error);
 215     assert(cl_error == CL_SUCCESS);
 216     ad->bShiftVecUploaded = false;
 217
 218     /* An element of the fshift device buffer has the same size as one element
 219        of the host side fshift buffer. */
 220     ad->fshift_elem_size = sizeof(*(((cl_nb_staging_t*)0)->fshift));
 221
 222     ad->fshift = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, SHIFTS * ad->fshift_elem_size, NULL, &cl_error);
 223     assert(cl_error == CL_SUCCESS);
 224     // TODO: handle errors, check clCreateBuffer flags
 225
 226     ad->e_lj = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, sizeof(float), NULL, &cl_error);
 227     assert(cl_error == CL_SUCCESS);
 228     // TODO: handle errors, check clCreateBuffer flags
 229
 230     ad->e_el = clCreateBuffer(runData->context, CL_MEM_READ_WRITE, sizeof(float), NULL, &cl_error);
 231     assert(cl_error == CL_SUCCESS);
 232     // TODO: handle errors, check clCreateBuffer flags
 233
 234     /* initialize to NULL pointers to data that is not allocated here and will
 235        need reallocation in nbnxn_gpu_init_atomdata */
 236     ad->xq = NULL;
 237     ad->f  = NULL;
 238
 239     /* size -1 indicates that the respective array hasn't been initialized yet */
 240     ad->natoms = -1;
 241     ad->nalloc = -1;
 242 }
 243
 244 /*! \brief Copies all parameters related to the cut-off from ic to nbp
 245  */
 246 static void set_cutoff_parameters(cl_nbparam_t              *nbp,
 247                                   const interaction_const_t *ic,
 248                                   const NbnxnListParameters *listParams)
 249 {
 250     nbp->ewald_beta        = ic->ewaldcoeff_q;
 251     nbp->sh_ewald          = ic->sh_ewald;
 252     nbp->epsfac            = ic->epsfac;
 253     nbp->two_k_rf          = 2.0 * ic->k_rf;
 254     nbp->c_rf              = ic->c_rf;
 255     nbp->rvdw_sq           = ic->rvdw * ic->rvdw;
 256     nbp->rcoulomb_sq       = ic->rcoulomb * ic->rcoulomb;
 257     nbp->rlistOuter_sq     = listParams->rlistOuter * listParams->rlistOuter;
 258     nbp->rlistInner_sq     = listParams->rlistInner * listParams->rlistInner;
 259     nbp->useDynamicPruning = listParams->useDynamicPruning;
 260
 261     nbp->sh_lj_ewald       = ic->sh_lj_ewald;
 262     nbp->ewaldcoeff_lj     = ic->ewaldcoeff_lj;
 263
 264     nbp->rvdw_switch       = ic->rvdw_switch;
 265     nbp->dispersion_shift  = ic->dispersion_shift;
 266     nbp->repulsion_shift   = ic->repulsion_shift;
 267     nbp->vdw_switch        = ic->vdw_switch;
 268 }
 269
 270 /*! \brief Returns the kinds of electrostatics and Vdw OpenCL
 271  *  kernels that will be used.
 272  *
 273  * Respectively, these values are from enum eelOcl and enum
 274  * evdwOcl. */
 275 static void
 276 map_interaction_types_to_gpu_kernel_flavors(const interaction_const_t *ic,
 277                                             int                        combRule,
 278                                             int                       *gpu_eeltype,
 279                                             int                       *gpu_vdwtype)
 280 {
 281     if (ic->vdwtype == evdwCUT)
 282     {
 283         switch (ic->vdw_modifier)
 284         {
 285             case eintmodNONE:
 286             case eintmodPOTSHIFT:
 287                 switch (combRule)
 288                 {
 289                     case ljcrNONE:
 290                         *gpu_vdwtype = evdwOclCUT;
 291                         break;
 292                     case ljcrGEOM:
 293                         *gpu_vdwtype = evdwOclCUTCOMBGEOM;
 294                         break;
 295                     case ljcrLB:
 296                         *gpu_vdwtype = evdwOclCUTCOMBLB;
 297                         break;
 298                     default:
 299                         gmx_incons("The requested LJ combination rule is not implemented in the OpenCL GPU accelerated kernels!");
 300                         break;
 301                 }
 302                 break;
 303             case eintmodFORCESWITCH:
 304                 *gpu_vdwtype = evdwOclFSWITCH;
 305                 break;
 306             case eintmodPOTSWITCH:
 307                 *gpu_vdwtype = evdwOclPSWITCH;
 308                 break;
 309             default:
 310                 gmx_incons("The requested VdW interaction modifier is not implemented in the GPU accelerated kernels!");
 311                 break;
 312         }
 313     }
 314     else if (ic->vdwtype == evdwPME)
 315     {
 316         if (ic->ljpme_comb_rule == ljcrGEOM)
 317         {
 318             *gpu_vdwtype = evdwOclEWALDGEOM;
 319         }
 320         else
 321         {
 322             *gpu_vdwtype = evdwOclEWALDLB;
 323         }
 324     }
 325     else
 326     {
 327         gmx_incons("The requested VdW type is not implemented in the GPU accelerated kernels!");
 328     }
 329
 330     if (ic->eeltype == eelCUT)
 331     {
 332         *gpu_eeltype = eelOclCUT;
 333     }
 334     else if (EEL_RF(ic->eeltype))
 335     {
 336         *gpu_eeltype = eelOclRF;
 337     }
 338     else if ((EEL_PME(ic->eeltype) || ic->eeltype == eelEWALD))
 339     {
 340         /* Initially rcoulomb == rvdw, so it's surely not twin cut-off. */
 341         *gpu_eeltype = nbnxn_gpu_pick_ewald_kernel_type(false);
 342     }
 343     else
 344     {
 345         /* Shouldn't happen, as this is checked when choosing Verlet-scheme */
 346         gmx_incons("The requested electrostatics type is not implemented in the GPU accelerated kernels!");
 347     }
 348 }
 349
 350 /*! \brief Initializes the nonbonded parameter data structure.
 351  */
 352 static void init_nbparam(cl_nbparam_t                    *nbp,
 353                          const interaction_const_t       *ic,
 354                          const NbnxnListParameters       *listParams,
 355                          const nbnxn_atomdata_t          *nbat,
 356                          const gmx_device_runtime_data_t *runData)
 357 {
 358     int         ntypes, nnbfp, nnbfp_comb;
 359     cl_int      cl_error;
 360
 361
 362     ntypes = nbat->ntype;
 363
 364     set_cutoff_parameters(nbp, ic, listParams);
 365
 366     map_interaction_types_to_gpu_kernel_flavors(ic,
 367                                                 nbat->comb_rule,
 368                                                 &(nbp->eeltype),
 369                                                 &(nbp->vdwtype));
 370
 371     if (ic->vdwtype == evdwPME)
 372     {
 373         if (ic->ljpme_comb_rule == ljcrGEOM)
 374         {
 375             assert(nbat->comb_rule == ljcrGEOM);
 376         }
 377         else
 378         {
 379             assert(nbat->comb_rule == ljcrLB);
 380         }
 381     }
 382     /* generate table for PME */
 383     nbp->coulomb_tab_climg2d = NULL;
 384     if (nbp->eeltype == eelOclEWALD_TAB || nbp->eeltype == eelOclEWALD_TAB_TWIN)
 385     {
 386         init_ewald_coulomb_force_table(ic, nbp, runData);
 387     }
 388     else
 389     // TODO: improvement needed.
 390     // The image2d is created here even if eeltype is not eelCuEWALD_TAB or eelCuEWALD_TAB_TWIN because the OpenCL kernels
 391     // don't accept NULL values for image2D parameters.
 392     {
 393         /* Switched from using textures to using buffers */
 394         // TODO: decide which alternative is most efficient - textures or buffers.
 395         /*
 396            cl_image_format array_format;
 397
 398            array_format.image_channel_data_type = CL_FLOAT;
 399            array_format.image_channel_order     = CL_R;
 400
 401            nbp->coulomb_tab_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE,
 402             &array_format, 1, 1, 0, NULL, &cl_error);
 403          */
 404
 405         nbp->coulomb_tab_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY, sizeof(cl_float), NULL, &cl_error);
 406         // TODO: handle errors
 407     }
 408
 409     nnbfp      = 2*ntypes*ntypes;
 410     nnbfp_comb = 2*ntypes;
 411
 412     {
 413         /* Switched from using textures to using buffers */
 414         // TODO: decide which alternative is most efficient - textures or buffers.
 415         /*
 416            cl_image_format array_format;
 417
 418            array_format.image_channel_data_type = CL_FLOAT;
 419            array_format.image_channel_order     = CL_R;
 420
 421            nbp->nbfp_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,
 422             &array_format, nnbfp, 1, 0, nbat->nbfp, &cl_error);
 423          */
 424
 425         nbp->nbfp_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, nnbfp*sizeof(cl_float), nbat->nbfp, &cl_error);
 426         assert(cl_error == CL_SUCCESS);
 427         // TODO: handle errors
 428
 429         if (ic->vdwtype == evdwPME)
 430         {
 431             /* Switched from using textures to using buffers */
 432             // TODO: decide which alternative is most efficient - textures or buffers.
 433             /*  nbp->nbfp_comb_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR,
 434                 &array_format, nnbfp_comb, 1, 0, nbat->nbfp_comb, &cl_error);*/
 435             nbp->nbfp_comb_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, nnbfp_comb*sizeof(cl_float), nbat->nbfp_comb, &cl_error);
 436
 437
 438             assert(cl_error == CL_SUCCESS);
 439             // TODO: handle errors
 440         }
 441         else
 442         {
 443             // TODO: improvement needed.
 444             // The image2d is created here even if vdwtype is not evdwPME because the OpenCL kernels
 445             // don't accept NULL values for image2D parameters.
 446             /* Switched from using textures to using buffers */
 447             // TODO: decide which alternative is most efficient - textures or buffers.
 448             /* nbp->nbfp_comb_climg2d = clCreateImage2D(runData->context, CL_MEM_READ_WRITE,
 449                 &array_format, 1, 1, 0, NULL, &cl_error);*/
 450             nbp->nbfp_comb_climg2d = clCreateBuffer(runData->context, CL_MEM_READ_ONLY, sizeof(cl_float), NULL, &cl_error);
 451
 452
 453             assert(cl_error == CL_SUCCESS);
 454             // TODO: handle errors
 455         }
 456     }
 457 }
 458
 459 //! This function is documented in the header file
 460 void nbnxn_gpu_pme_loadbal_update_param(const nonbonded_verlet_t    *nbv,
 461                                         const interaction_const_t   *ic,
 462                                         const NbnxnListParameters   *listParams)
 463 {
 464     if (!nbv || nbv->grp[0].kernel_type != nbnxnk8x8x8_GPU)
 465     {
 466         return;
 467     }
 468     gmx_nbnxn_ocl_t    *nb  = nbv->gpu_nbv;
 469     cl_nbparam_t       *nbp = nb->nbparam;
 470
 471     set_cutoff_parameters(nbp, ic, listParams);
 472
 473     nbp->eeltype = nbnxn_gpu_pick_ewald_kernel_type(ic->rcoulomb != ic->rvdw);
 474
 475     init_ewald_coulomb_force_table(ic, nb->nbparam, nb->dev_rundata);
 476 }
 477
 478 /*! \brief Initializes the pair list data structure.
 479  */
 480 static void init_plist(cl_plist_t *pl)
 481 {
 482     /* initialize to NULL pointers to data that is not allocated here and will
 483        need reallocation in nbnxn_gpu_init_pairlist */
 484     pl->sci     = NULL;
 485     pl->cj4     = NULL;
 486     pl->imask   = NULL;
 487     pl->excl    = NULL;
 488
 489     /* size -1 indicates that the respective array hasn't been initialized yet */
 490     pl->na_c           = -1;
 491     pl->nsci           = -1;
 492     pl->sci_nalloc     = -1;
 493     pl->ncj4           = -1;
 494     pl->cj4_nalloc     = -1;
 495     pl->nimask         = -1;
 496     pl->imask_nalloc   = -1;
 497     pl->nexcl          = -1;
 498     pl->excl_nalloc    = -1;
 499     pl->haveFreshList  = false;
 500 }
 501
 502 /*! \brief Initializes the timer data structure.
 503  */
 504 static void init_timers(cl_timers_t *t,
 505                         bool         bUseTwoStreams)
 506 {
 507     for (int i = 0; i <= (bUseTwoStreams ? 1 : 0); i++)
 508     {
 509         t->didPairlistH2D[i]  = false;
 510         t->didPrune[i]        = false;
 511         t->didRollingPrune[i] = false;
 512     }
 513 }
 514
 515 /*! \brief Initializes the timings data structure.
 516  */
 517 static void init_timings(gmx_wallclock_gpu_nbnxn_t *t)
 518 {
 519     int i, j;
 520
 521     t->nb_h2d_t = 0.0;
 522     t->nb_d2h_t = 0.0;
 523     t->nb_c     = 0;
 524     t->pl_h2d_t = 0.0;
 525     t->pl_h2d_c = 0;
 526     for (i = 0; i < 2; i++)
 527     {
 528         for (j = 0; j < 2; j++)
 529         {
 530             t->ktime[i][j].t = 0.0;
 531             t->ktime[i][j].c = 0;
 532         }
 533     }
 534
 535     t->pruneTime.c        = 0;
 536     t->pruneTime.t        = 0.0;
 537     t->dynamicPruneTime.c = 0;
 538     t->dynamicPruneTime.t = 0.0;
 539 }
 540
 541 /*! \brief Creates context for OpenCL GPU given by \p mygpu
 542  *
 543  * A fatal error results if creation fails.
 544  *
 545  * \param[inout] runtimeData runtime data including program and context
 546  * \param[in]    devInfo     device info struct
 547  * \param[in]    rank        MPI rank (for error reporting)
 548  */
 549 static void
 550 nbnxn_gpu_create_context(gmx_device_runtime_data_t *runtimeData,
 551                          const gmx_device_info_t   *devInfo,
 552                          int                        rank)
 553 {
 554     cl_context_properties     context_properties[3];
 555     cl_platform_id            platform_id;
 556     cl_device_id              device_id;
 557     cl_context                context;
 558     cl_int                    cl_error;
 559
 560     assert(runtimeData != NULL);
 561     assert(devInfo != NULL);
 562
 563     platform_id      = devInfo->ocl_gpu_id.ocl_platform_id;
 564     device_id        = devInfo->ocl_gpu_id.ocl_device_id;
 565
 566     context_properties[0] = CL_CONTEXT_PLATFORM;
 567     context_properties[1] = (cl_context_properties) platform_id;
 568     context_properties[2] = 0; /* Terminates the list of properties */
 569
 570     context = clCreateContext(context_properties, 1, &device_id, NULL, NULL, &cl_error);
 571     if (CL_SUCCESS != cl_error)
 572     {
 573         gmx_fatal(FARGS, "On rank %d failed to create context for GPU #%s:\n OpenCL error %d: %s",
 574                   rank,
 575                   devInfo->device_name,
 576                   cl_error, ocl_get_error_string(cl_error).c_str());
 577         return;
 578     }
 579
 580     runtimeData->context = context;
 581 }
 582
 583 /*! \brief Initializes the OpenCL kernel pointers of the nbnxn_ocl_ptr_t input data structure. */
 584 static cl_kernel nbnxn_gpu_create_kernel(gmx_nbnxn_ocl_t *nb,
 585                                          const char      *kernel_name)
 586 {
 587     cl_kernel kernel;
 588     cl_int    cl_error;
 589
 590     kernel = clCreateKernel(nb->dev_rundata->program, kernel_name, &cl_error);
 591     if (CL_SUCCESS != cl_error)
 592     {
 593         gmx_fatal(FARGS, "Failed to create kernel '%s' for GPU #%s: OpenCL error %d",
 594                   kernel_name,
 595                   nb->dev_info->device_name,
 596                   cl_error);
 597     }
 598
 599     return kernel;
 600 }
 601
 602 /*! \brief Clears nonbonded shift force output array and energy outputs on the GPU.
 603  */
 604 static void
 605 nbnxn_ocl_clear_e_fshift(gmx_nbnxn_ocl_t *nb)
 606 {
 607
 608     cl_int               cl_error;
 609     cl_atomdata_t *      adat     = nb->atdat;
 610     cl_command_queue     ls       = nb->stream[eintLocal];
 611
 612     size_t               local_work_size[3]   = {1, 1, 1};
 613     size_t               global_work_size[3]  = {1, 1, 1};
 614
 615     cl_int               shifts   = SHIFTS*3;
 616
 617     cl_int               arg_no;
 618
 619     cl_kernel            zero_e_fshift = nb->kernel_zero_e_fshift;
 620
 621     local_work_size[0]   = 64;
 622     // Round the total number of threads up from the array size
 623     global_work_size[0]  = ((shifts + local_work_size[0] - 1)/local_work_size[0])*local_work_size[0];
 624
 625     arg_no    = 0;
 626     cl_error  = clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_mem), &(adat->fshift));
 627     cl_error |= clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_mem), &(adat->e_lj));
 628     cl_error |= clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_mem), &(adat->e_el));
 629     cl_error |= clSetKernelArg(zero_e_fshift, arg_no++, sizeof(cl_uint), &shifts);
 630     GMX_ASSERT(cl_error == CL_SUCCESS, ocl_get_error_string(cl_error).c_str());
 631
 632     cl_error = clEnqueueNDRangeKernel(ls, zero_e_fshift, 3, NULL, global_work_size, local_work_size, 0, NULL, NULL);
 633     GMX_ASSERT(cl_error == CL_SUCCESS, ocl_get_error_string(cl_error).c_str());
 634 }
 635
 636 /*! \brief Initializes the OpenCL kernel pointers of the nbnxn_ocl_ptr_t input data structure. */
 637 static void nbnxn_gpu_init_kernels(gmx_nbnxn_ocl_t *nb)
 638 {
 639     /* Init to 0 main kernel arrays */
 640     /* They will be later on initialized in select_nbnxn_kernel */
 641     // TODO: consider always creating all variants of the kernels here so that there is no
 642     // need for late call to clCreateKernel -- if that gives any advantage?
 643     memset(nb->kernel_ener_noprune_ptr, 0, sizeof(nb->kernel_ener_noprune_ptr));
 644     memset(nb->kernel_ener_prune_ptr, 0, sizeof(nb->kernel_ener_prune_ptr));
 645     memset(nb->kernel_noener_noprune_ptr, 0, sizeof(nb->kernel_noener_noprune_ptr));
 646     memset(nb->kernel_noener_prune_ptr, 0, sizeof(nb->kernel_noener_prune_ptr));
 647
 648     /* Init pruning kernels
 649      *
 650      * TODO: we could avoid creating kernels if dynamic pruning is turned off,
 651      * but ATM that depends on force flags not passed into the initialization.
 652      */
 653     nb->kernel_pruneonly[epruneFirst]   = nbnxn_gpu_create_kernel(nb, "nbnxn_kernel_prune_opencl");
 654     nb->kernel_pruneonly[epruneRolling] = nbnxn_gpu_create_kernel(nb, "nbnxn_kernel_prune_rolling_opencl");
 655
 656     /* Init auxiliary kernels */
 657     nb->kernel_memset_f      = nbnxn_gpu_create_kernel(nb, "memset_f");
 658     nb->kernel_memset_f2     = nbnxn_gpu_create_kernel(nb, "memset_f2");
 659     nb->kernel_memset_f3     = nbnxn_gpu_create_kernel(nb, "memset_f3");
 660     nb->kernel_zero_e_fshift = nbnxn_gpu_create_kernel(nb, "zero_e_fshift");
 661 }
 662
 663 /*! \brief Initializes simulation constant data.
 664  *
 665  *  Initializes members of the atomdata and nbparam structs and
 666  *  clears e/fshift output buffers.
 667  */
 668 static void nbnxn_ocl_init_const(gmx_nbnxn_ocl_t                *nb,
 669                                  const interaction_const_t      *ic,
 670                                  const NbnxnListParameters      *listParams,
 671                                  const nbnxn_atomdata_t         *nbat)
 672 {
 673     init_atomdata_first(nb->atdat, nbat->ntype, nb->dev_rundata);
 674     init_nbparam(nb->nbparam, ic, listParams, nbat, nb->dev_rundata);
 675 }
 676
 677
 678 //! This function is documented in the header file
 679 void nbnxn_gpu_init(gmx_nbnxn_ocl_t          **p_nb,
 680                     const gmx_device_info_t   *deviceInfo,
 681                     const interaction_const_t *ic,
 682                     const NbnxnListParameters *listParams,
 683                     const nbnxn_atomdata_t    *nbat,
 684                     int                        rank,
 685                     gmx_bool                   bLocalAndNonlocal)
 686 {
 687     gmx_nbnxn_ocl_t            *nb;
 688     cl_int                      cl_error;
 689     cl_command_queue_properties queue_properties;
 690
 691     assert(ic);
 692
 693     if (p_nb == NULL)
 694     {
 695         return;
 696     }
 697
 698     snew(nb, 1);
 699     snew(nb->atdat, 1);
 700     snew(nb->nbparam, 1);
 701     snew(nb->plist[eintLocal], 1);
 702     if (bLocalAndNonlocal)
 703     {
 704         snew(nb->plist[eintNonlocal], 1);
 705     }
 706
 707     nb->bUseTwoStreams = bLocalAndNonlocal;
 708
 709     nb->timers = new cl_timers_t();
 710     snew(nb->timings, 1);
 711
 712     /* set device info, just point it to the right GPU among the detected ones */
 713     nb->dev_info = deviceInfo;
 714     snew(nb->dev_rundata, 1);
 715
 716     /* init to NULL the debug buffer */
 717     nb->debug_buffer = NULL;
 718
 719     /* init nbst */
 720     ocl_pmalloc((void**)&nb->nbst.e_lj, sizeof(*nb->nbst.e_lj));
 721     ocl_pmalloc((void**)&nb->nbst.e_el, sizeof(*nb->nbst.e_el));
 722     ocl_pmalloc((void**)&nb->nbst.fshift, SHIFTS * sizeof(*nb->nbst.fshift));
 723
 724     init_plist(nb->plist[eintLocal]);
 725
 726     /* OpenCL timing disabled if GMX_DISABLE_GPU_TIMING is defined. */
 727     nb->bDoTime = (getenv("GMX_DISABLE_GPU_TIMING") == NULL);
 728
 729     /* Create queues only after bDoTime has been initialized */
 730     if (nb->bDoTime)
 731     {
 732         queue_properties = CL_QUEUE_PROFILING_ENABLE;
 733     }
 734     else
 735     {
 736         queue_properties = 0;
 737     }
 738
 739     nbnxn_gpu_create_context(nb->dev_rundata, nb->dev_info, rank);
 740
 741     /* local/non-local GPU streams */
 742     nb->stream[eintLocal] = clCreateCommandQueue(nb->dev_rundata->context, nb->dev_info->ocl_gpu_id.ocl_device_id, queue_properties, &cl_error);
 743     if (CL_SUCCESS != cl_error)
 744     {
 745         gmx_fatal(FARGS, "On rank %d failed to create context for GPU #%s: OpenCL error %d",
 746                   rank,
 747                   nb->dev_info->device_name,
 748                   cl_error);
 749         return;
 750     }
 751
 752     if (nb->bUseTwoStreams)
 753     {
 754         init_plist(nb->plist[eintNonlocal]);
 755
 756         nb->stream[eintNonlocal] = clCreateCommandQueue(nb->dev_rundata->context, nb->dev_info->ocl_gpu_id.ocl_device_id, queue_properties, &cl_error);
 757         if (CL_SUCCESS != cl_error)
 758         {
 759             gmx_fatal(FARGS, "On rank %d failed to create context for GPU #%s: OpenCL error %d",
 760                       rank,
 761                       nb->dev_info->device_name,
 762                       cl_error);
 763             return;
 764         }
 765     }
 766
 767     if (nb->bDoTime)
 768     {
 769         init_timers(nb->timers, nb->bUseTwoStreams);
 770         init_timings(nb->timings);
 771     }
 772
 773     nbnxn_ocl_init_const(nb, ic, listParams, nbat);
 774
 775     /* Enable LJ param manual prefetch for AMD or if we request through env. var.
 776      * TODO: decide about NVIDIA
 777      */
 778     nb->bPrefetchLjParam =
 779         (getenv("GMX_OCL_DISABLE_I_PREFETCH") == NULL) &&
 780         ((nb->dev_info->vendor_e == OCL_VENDOR_AMD) || (getenv("GMX_OCL_ENABLE_I_PREFETCH") != NULL));
 781
 782     /* NOTE: in CUDA we pick L1 cache configuration for the nbnxn kernels here,
 783      * but sadly this is not supported in OpenCL (yet?). Consider adding it if
 784      * it becomes supported.
 785      */
 786     nbnxn_gpu_compile_kernels(nb);
 787     nbnxn_gpu_init_kernels(nb);
 788
 789     /* clear energy and shift force outputs */
 790     nbnxn_ocl_clear_e_fshift(nb);
 791
 792     *p_nb = nb;
 793
 794     if (debug)
 795     {
 796         fprintf(debug, "Initialized OpenCL data structures.\n");
 797     }
 798 }
 799
 800 /*! \brief Clears the first natoms_clear elements of the GPU nonbonded force output array.
 801  */
 802 static void nbnxn_ocl_clear_f(gmx_nbnxn_ocl_t *nb, int natoms_clear)
 803 {
 804     if (natoms_clear == 0)
 805     {
 806         return;
 807     }
 808
 809     cl_int               cl_error;
 810     cl_atomdata_t *      adat     = nb->atdat;
 811     cl_command_queue     ls       = nb->stream[eintLocal];
 812     cl_float             value    = 0.0f;
 813
 814     size_t               local_work_size[3]  = {1, 1, 1};
 815     size_t               global_work_size[3] = {1, 1, 1};
 816
 817     cl_int               arg_no;
 818
 819     cl_kernel            memset_f = nb->kernel_memset_f;
 820
 821     cl_uint              natoms_flat = natoms_clear * (sizeof(rvec)/sizeof(real));
 822
 823     local_work_size[0]  = 64;
 824     // Round the total number of threads up from the array size
 825     global_work_size[0] = ((natoms_flat + local_work_size[0] - 1)/local_work_size[0])*local_work_size[0];
 826
 827
 828     arg_no    = 0;
 829     cl_error  = clSetKernelArg(memset_f, arg_no++, sizeof(cl_mem), &(adat->f));
 830     cl_error |= clSetKernelArg(memset_f, arg_no++, sizeof(cl_float), &value);
 831     cl_error |= clSetKernelArg(memset_f, arg_no++, sizeof(cl_uint), &natoms_flat);
 832     assert(cl_error == CL_SUCCESS);
 833
 834     cl_error = clEnqueueNDRangeKernel(ls, memset_f, 3, NULL, global_work_size, local_work_size, 0, NULL, NULL);
 835     assert(cl_error == CL_SUCCESS);
 836 }
 837
 838 //! This function is documented in the header file
 839 void
 840 nbnxn_gpu_clear_outputs(gmx_nbnxn_ocl_t   *nb,
 841                         int                flags)
 842 {
 843     nbnxn_ocl_clear_f(nb, nb->atdat->natoms);
 844     /* clear shift force array and energies if the outputs were
 845        used in the current step */
 846     if (flags & GMX_FORCE_VIRIAL)
 847     {
 848         nbnxn_ocl_clear_e_fshift(nb);
 849     }
 850
 851     /* kick off buffer clearing kernel to ensure concurrency with constraints/update */
 852     cl_int gmx_unused cl_error;
 853     cl_error = clFlush(nb->stream[eintLocal]);
 854     assert(CL_SUCCESS == cl_error);
 855 }
 856
 857 //! This function is documented in the header file
 858 void nbnxn_gpu_init_pairlist(gmx_nbnxn_ocl_t        *nb,
 859                              const nbnxn_pairlist_t *h_plist,
 860                              int                     iloc)
 861 {
 862     if (canSkipWork(nb, iloc))
 863     {
 864         return;
 865     }
 866
 867     char             sbuf[STRLEN];
 868     bool             bDoTime    = nb->bDoTime;
 869     cl_command_queue stream     = nb->stream[iloc];
 870     cl_plist_t      *d_plist    = nb->plist[iloc];
 871
 872     if (d_plist->na_c < 0)
 873     {
 874         d_plist->na_c = h_plist->na_ci;
 875     }
 876     else
 877     {
 878         if (d_plist->na_c != h_plist->na_ci)
 879         {
 880             sprintf(sbuf, "In cu_init_plist: the #atoms per cell has changed (from %d to %d)",
 881                     d_plist->na_c, h_plist->na_ci);
 882             gmx_incons(sbuf);
 883         }
 884     }
 885
 886     if (bDoTime)
 887     {
 888         nb->timers->pl_h2d[iloc].openTimingRegion(stream);
 889         nb->timers->didPairlistH2D[iloc] = true;
 890     }
 891
 892     ocl_realloc_buffered(&d_plist->sci, h_plist->sci, sizeof(nbnxn_sci_t),
 893                          &d_plist->nsci, &d_plist->sci_nalloc,
 894                          h_plist->nsci,
 895                          nb->dev_rundata->context,
 896                          stream, true, bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
 897
 898     ocl_realloc_buffered(&d_plist->cj4, h_plist->cj4, sizeof(nbnxn_cj4_t),
 899                          &d_plist->ncj4, &d_plist->cj4_nalloc,
 900                          h_plist->ncj4,
 901                          nb->dev_rundata->context,
 902                          stream, true, bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
 903
 904     /* this call only allocates space on the device (no data is transferred) - no timing as well! */
 905     ocl_realloc_buffered(&d_plist->imask, NULL, sizeof(unsigned int),
 906                          &d_plist->nimask, &d_plist->imask_nalloc,
 907                          h_plist->ncj4*c_nbnxnGpuClusterpairSplit,
 908                          nb->dev_rundata->context,
 909                          stream, true);
 910
 911     ocl_realloc_buffered(&d_plist->excl, h_plist->excl, sizeof(nbnxn_excl_t),
 912                          &d_plist->nexcl, &d_plist->excl_nalloc,
 913                          h_plist->nexcl,
 914                          nb->dev_rundata->context,
 915                          stream, true, bDoTime ? nb->timers->pl_h2d[iloc].fetchNextEvent() : nullptr);
 916
 917     if (bDoTime)
 918     {
 919         nb->timers->pl_h2d[iloc].closeTimingRegion(stream);
 920     }
 921
 922     /* need to prune the pair list during the next step */
 923     d_plist->haveFreshList = true;
 924 }
 925
 926 //! This function is documented in the header file
 927 void nbnxn_gpu_upload_shiftvec(gmx_nbnxn_ocl_t        *nb,
 928                                const nbnxn_atomdata_t *nbatom)
 929 {
 930     cl_atomdata_t   *adat  = nb->atdat;
 931     cl_command_queue ls    = nb->stream[eintLocal];
 932
 933     /* only if we have a dynamic box */
 934     if (nbatom->bDynamicBox || !adat->bShiftVecUploaded)
 935     {
 936         ocl_copy_H2D_async(adat->shift_vec, nbatom->shift_vec, 0,
 937                            SHIFTS * adat->shift_vec_elem_size, ls, NULL);
 938         adat->bShiftVecUploaded = true;
 939     }
 940 }
 941
 942 //! This function is documented in the header file
 943 void nbnxn_gpu_init_atomdata(gmx_nbnxn_ocl_t               *nb,
 944                              const nbnxn_atomdata_t        *nbat)
 945 {
 946     cl_int           cl_error;
 947     int              nalloc, natoms;
 948     bool             realloced;
 949     bool             bDoTime = nb->bDoTime;
 950     cl_timers_t     *timers  = nb->timers;
 951     cl_atomdata_t   *d_atdat = nb->atdat;
 952     cl_command_queue ls      = nb->stream[eintLocal];
 953
 954     natoms    = nbat->natoms;
 955     realloced = false;
 956
 957     if (bDoTime)
 958     {
 959         /* time async copy */
 960         timers->atdat.openTimingRegion(ls);
 961     }
 962
 963     /* need to reallocate if we have to copy more atoms than the amount of space
 964        available and only allocate if we haven't initialized yet, i.e d_atdat->natoms == -1 */
 965     if (natoms > d_atdat->nalloc)
 966     {
 967         nalloc = over_alloc_small(natoms);
 968
 969         /* free up first if the arrays have already been initialized */
 970         if (d_atdat->nalloc != -1)
 971         {
 972             freeDeviceBuffer(&d_atdat->f);
 973             freeDeviceBuffer(&d_atdat->xq);
 974             freeDeviceBuffer(&d_atdat->lj_comb);
 975             freeDeviceBuffer(&d_atdat->atom_types);
 976         }
 977
 978         d_atdat->f_elem_size = sizeof(rvec);
 979
 980         // TODO: handle errors, check clCreateBuffer flags
 981         d_atdat->f = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * d_atdat->f_elem_size, NULL, &cl_error);
 982         assert(CL_SUCCESS == cl_error);
 983
 984         // TODO: change the flag to read-only
 985         d_atdat->xq = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * sizeof(cl_float4), NULL, &cl_error);
 986         assert(CL_SUCCESS == cl_error);
 987         // TODO: handle errors, check clCreateBuffer flags
 988
 989         if (useLjCombRule(nb->nbparam->vdwtype))
 990         {
 991             // TODO: change the flag to read-only
 992             d_atdat->lj_comb = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * sizeof(cl_float2), NULL, &cl_error);
 993             assert(CL_SUCCESS == cl_error);
 994             // TODO: handle errors, check clCreateBuffer flags
 995         }
 996         else
 997         {
 998             // TODO: change the flag to read-only
 999             d_atdat->atom_types = clCreateBuffer(nb->dev_rundata->context, CL_MEM_READ_WRITE, nalloc * sizeof(int), NULL, &cl_error);
1000             assert(CL_SUCCESS == cl_error);
1001             // TODO: handle errors, check clCreateBuffer flags
1002         }
1003
1004         d_atdat->nalloc = nalloc;
1005         realloced       = true;
1006     }
1007
1008     d_atdat->natoms       = natoms;
1009     d_atdat->natoms_local = nbat->natoms_local;
1010
1011     /* need to clear GPU f output if realloc happened */
1012     if (realloced)
1013     {
1014         nbnxn_ocl_clear_f(nb, nalloc);
1015     }
1016
1017     if (useLjCombRule(nb->nbparam->vdwtype))
1018     {
1019         ocl_copy_H2D_async(d_atdat->lj_comb, nbat->lj_comb, 0,
1020                            natoms*sizeof(cl_float2), ls, bDoTime ? timers->atdat.fetchNextEvent() : nullptr);
1021     }
1022     else
1023     {
1024         ocl_copy_H2D_async(d_atdat->atom_types, nbat->type, 0,
1025                            natoms*sizeof(int), ls, bDoTime ? timers->atdat.fetchNextEvent() : nullptr);
1026
1027     }
1028
1029     if (bDoTime)
1030     {
1031         timers->atdat.closeTimingRegion(ls);
1032     }
1033
1034     /* kick off the tasks enqueued above to ensure concurrency with the search */
1035     cl_error = clFlush(ls);
1036     assert(CL_SUCCESS == cl_error);
1037 }
1038
1039 /*! \brief Releases an OpenCL kernel pointer */
1040 static void free_kernel(cl_kernel *kernel_ptr)
1041 {
1042     cl_int gmx_unused cl_error;
1043
1044     assert(NULL != kernel_ptr);
1045
1046     if (*kernel_ptr)
1047     {
1048         cl_error = clReleaseKernel(*kernel_ptr);
1049         assert(cl_error == CL_SUCCESS);
1050
1051         *kernel_ptr = NULL;
1052     }
1053 }
1054
1055 /*! \brief Releases a list of OpenCL kernel pointers */
1056 static void free_kernels(cl_kernel *kernels, int count)
1057 {
1058     int i;
1059
1060     for (i = 0; i < count; i++)
1061     {
1062         free_kernel(kernels + i);
1063     }
1064 }
1065
1066 /*! \brief Free the OpenCL runtime data (context and program).
1067  *
1068  *  The function releases the OpenCL context and program assuciated with the
1069  *  device that the calling PP rank is running on.
1070  *
1071  *  \param runData [in]  porinter to the structure with runtime data.
1072  */
1073 static void free_gpu_device_runtime_data(gmx_device_runtime_data_t *runData)
1074 {
1075     if (runData == NULL)
1076     {
1077         return;
1078     }
1079
1080     cl_int gmx_unused cl_error;
1081
1082     if (runData->context)
1083     {
1084         cl_error         = clReleaseContext(runData->context);
1085         runData->context = NULL;
1086         assert(CL_SUCCESS == cl_error);
1087     }
1088
1089     if (runData->program)
1090     {
1091         cl_error         = clReleaseProgram(runData->program);
1092         runData->program = NULL;
1093         assert(CL_SUCCESS == cl_error);
1094     }
1095
1096 }
1097
1098 //! This function is documented in the header file
1099 void nbnxn_gpu_free(gmx_nbnxn_ocl_t *nb)
1100 {
1101     if (nb == NULL)
1102     {
1103         return;
1104     }
1105
1106     /* Free kernels */
1107     int kernel_count = sizeof(nb->kernel_ener_noprune_ptr) / sizeof(nb->kernel_ener_noprune_ptr[0][0]);
1108     free_kernels((cl_kernel*)nb->kernel_ener_noprune_ptr, kernel_count);
1109
1110     kernel_count = sizeof(nb->kernel_ener_prune_ptr) / sizeof(nb->kernel_ener_prune_ptr[0][0]);
1111     free_kernels((cl_kernel*)nb->kernel_ener_prune_ptr, kernel_count);
1112
1113     kernel_count = sizeof(nb->kernel_noener_noprune_ptr) / sizeof(nb->kernel_noener_noprune_ptr[0][0]);
1114     free_kernels((cl_kernel*)nb->kernel_noener_noprune_ptr, kernel_count);
1115
1116     kernel_count = sizeof(nb->kernel_noener_prune_ptr) / sizeof(nb->kernel_noener_prune_ptr[0][0]);
1117     free_kernels((cl_kernel*)nb->kernel_noener_prune_ptr, kernel_count);
1118
1119     free_kernel(&(nb->kernel_memset_f));
1120     free_kernel(&(nb->kernel_memset_f2));
1121     free_kernel(&(nb->kernel_memset_f3));
1122     free_kernel(&(nb->kernel_zero_e_fshift));
1123
1124     /* Free atdat */
1125     freeDeviceBuffer(&(nb->atdat->xq));
1126     freeDeviceBuffer(&(nb->atdat->f));
1127     freeDeviceBuffer(&(nb->atdat->e_lj));
1128     freeDeviceBuffer(&(nb->atdat->e_el));
1129     freeDeviceBuffer(&(nb->atdat->fshift));
1130     freeDeviceBuffer(&(nb->atdat->lj_comb));
1131     freeDeviceBuffer(&(nb->atdat->atom_types));
1132     freeDeviceBuffer(&(nb->atdat->shift_vec));
1133     sfree(nb->atdat);
1134
1135     /* Free nbparam */
1136     freeDeviceBuffer(&(nb->nbparam->nbfp_climg2d));
1137     freeDeviceBuffer(&(nb->nbparam->nbfp_comb_climg2d));
1138     freeDeviceBuffer(&(nb->nbparam->coulomb_tab_climg2d));
1139     sfree(nb->nbparam);
1140
1141     /* Free plist */
1142     auto *plist = nb->plist[eintLocal];
1143     freeDeviceBuffer(&plist->sci);
1144     freeDeviceBuffer(&plist->cj4);
1145     freeDeviceBuffer(&plist->imask);
1146     freeDeviceBuffer(&plist->excl);
1147     sfree(plist);
1148     if (nb->bUseTwoStreams)
1149     {
1150         auto *plist_nl = nb->plist[eintNonlocal];
1151         freeDeviceBuffer(&plist_nl->sci);
1152         freeDeviceBuffer(&plist_nl->cj4);
1153         freeDeviceBuffer(&plist_nl->imask);
1154         freeDeviceBuffer(&plist_nl->excl);
1155         sfree(plist_nl);
1156     }
1157
1158     /* Free nbst */
1159     ocl_pfree(nb->nbst.e_lj);
1160     nb->nbst.e_lj = NULL;
1161
1162     ocl_pfree(nb->nbst.e_el);
1163     nb->nbst.e_el = NULL;
1164
1165     ocl_pfree(nb->nbst.fshift);
1166     nb->nbst.fshift = NULL;
1167
1168     /* Free debug buffer */
1169     freeDeviceBuffer(&nb->debug_buffer);
1170
1171     /* Free command queues */
1172     clReleaseCommandQueue(nb->stream[eintLocal]);
1173     nb->stream[eintLocal] = NULL;
1174     if (nb->bUseTwoStreams)
1175     {
1176         clReleaseCommandQueue(nb->stream[eintNonlocal]);
1177         nb->stream[eintNonlocal] = NULL;
1178     }
1179     /* Free other events */
1180     if (nb->nonlocal_done)
1181     {
1182         clReleaseEvent(nb->nonlocal_done);
1183         nb->nonlocal_done = NULL;
1184     }
1185     if (nb->misc_ops_and_local_H2D_done)
1186     {
1187         clReleaseEvent(nb->misc_ops_and_local_H2D_done);
1188         nb->misc_ops_and_local_H2D_done = NULL;
1189     }
1190
1191     free_gpu_device_runtime_data(nb->dev_rundata);
1192     sfree(nb->dev_rundata);
1193
1194     /* Free timers and timings */
1195     delete nb->timers;
1196     sfree(nb->timings);
1197     sfree(nb);
1198
1199     if (debug)
1200     {
1201         fprintf(debug, "Cleaned up OpenCL data structures.\n");
1202     }
1203 }
1204
1205 //! This function is documented in the header file
1206 gmx_wallclock_gpu_nbnxn_t *nbnxn_gpu_get_timings(gmx_nbnxn_ocl_t *nb)
1207 {
1208     return (nb != nullptr && nb->bDoTime) ? nb->timings : nullptr;
1209 }
1210
1211 //! This function is documented in the header file
1212 void nbnxn_gpu_reset_timings(nonbonded_verlet_t* nbv)
1213 {
1214     if (nbv->gpu_nbv && nbv->gpu_nbv->bDoTime)
1215     {
1216         init_timings(nbv->gpu_nbv->timings);
1217     }
1218 }
1219
1220 //! This function is documented in the header file
1221 int nbnxn_gpu_min_ci_balanced(gmx_nbnxn_ocl_t *nb)
1222 {
1223     return nb != NULL ?
1224            gpu_min_ci_balanced_factor * nb->dev_info->compute_units : 0;
1225 }
1226
1227 //! This function is documented in the header file
1228 gmx_bool nbnxn_gpu_is_kernel_ewald_analytical(const gmx_nbnxn_ocl_t *nb)
1229 {
1230     return ((nb->nbparam->eeltype == eelOclEWALD_ANA) ||
1231             (nb->nbparam->eeltype == eelOclEWALD_ANA_TWIN));
1232 }