src/gromacs/mdlib/nbnxn_search.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2012, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  */
  32
  33 #ifdef HAVE_CONFIG_H
  34 #include <config.h>
  35 #endif
  36
  37 #include <math.h>
  38 #include <string.h>
  39 #include "sysstuff.h"
  40 #include "smalloc.h"
  41 #include "macros.h"
  42 #include "maths.h"
  43 #include "vec.h"
  44 #include "pbc.h"
  45 #include "nbnxn_consts.h"
  46 /* nbnxn_internal.h included gromacs/simd/macros.h */
  47 #include "nbnxn_internal.h"
  48 #ifdef GMX_NBNXN_SIMD
  49 #include "gromacs/simd/vector_operations.h"
  50 #endif
  51 #include "nbnxn_atomdata.h"
  52 #include "nbnxn_search.h"
  53 #include "gmx_omp_nthreads.h"
  54 #include "nrnb.h"
  55
  56 #include "gromacs/fileio/gmxfio.h"
  57
  58 #ifdef NBNXN_SEARCH_BB_SIMD4
  59 /* We use 4-wide SIMD for bounding box calculations */
  60
  61 #ifndef GMX_DOUBLE
  62 /* Single precision BBs + coordinates, we can also load coordinates with SIMD */
  63 #define NBNXN_SEARCH_SIMD4_FLOAT_X_BB
  64 #endif
  65
  66 #if defined NBNXN_SEARCH_SIMD4_FLOAT_X_BB && (GPU_NSUBCELL == 4 || GPU_NSUBCELL == 8)
  67 /* Store bounding boxes with x, y and z coordinates in packs of 4 */
  68 #define NBNXN_PBB_SIMD4
  69 #endif
  70
  71 /* The packed bounding box coordinate stride is always set to 4.
  72  * With AVX we could use 8, but that turns out not to be faster.
  73  */
  74 #define STRIDE_PBB        4
  75 #define STRIDE_PBB_2LOG   2
  76
  77 #endif /* NBNXN_SEARCH_BB_SIMD4 */
  78
  79 #ifdef GMX_NBNXN_SIMD
  80
  81 /* The functions below are macros as they are performance sensitive */
  82
  83 /* 4x4 list, pack=4: no complex conversion required */
  84 /* i-cluster to j-cluster conversion */
  85 #define CI_TO_CJ_J4(ci)   (ci)
  86 /* cluster index to coordinate array index conversion */
  87 #define X_IND_CI_J4(ci)  ((ci)*STRIDE_P4)
  88 #define X_IND_CJ_J4(cj)  ((cj)*STRIDE_P4)
  89
  90 /* 4x2 list, pack=4: j-cluster size is half the packing width */
  91 /* i-cluster to j-cluster conversion */
  92 #define CI_TO_CJ_J2(ci)  ((ci)<<1)
  93 /* cluster index to coordinate array index conversion */
  94 #define X_IND_CI_J2(ci)  ((ci)*STRIDE_P4)
  95 #define X_IND_CJ_J2(cj)  (((cj)>>1)*STRIDE_P4 + ((cj) & 1)*(PACK_X4>>1))
  96
  97 /* 4x8 list, pack=8: i-cluster size is half the packing width */
  98 /* i-cluster to j-cluster conversion */
  99 #define CI_TO_CJ_J8(ci)  ((ci)>>1)
 100 /* cluster index to coordinate array index conversion */
 101 #define X_IND_CI_J8(ci)  (((ci)>>1)*STRIDE_P8 + ((ci) & 1)*(PACK_X8>>1))
 102 #define X_IND_CJ_J8(cj)  ((cj)*STRIDE_P8)
 103
 104 /* The j-cluster size is matched to the SIMD width */
 105 #if GMX_SIMD_WIDTH_HERE == 2
 106 #define CI_TO_CJ_SIMD_4XN(ci)  CI_TO_CJ_J2(ci)
 107 #define X_IND_CI_SIMD_4XN(ci)  X_IND_CI_J2(ci)
 108 #define X_IND_CJ_SIMD_4XN(cj)  X_IND_CJ_J2(cj)
 109 #else
 110 #if GMX_SIMD_WIDTH_HERE == 4
 111 #define CI_TO_CJ_SIMD_4XN(ci)  CI_TO_CJ_J4(ci)
 112 #define X_IND_CI_SIMD_4XN(ci)  X_IND_CI_J4(ci)
 113 #define X_IND_CJ_SIMD_4XN(cj)  X_IND_CJ_J4(cj)
 114 #else
 115 #if GMX_SIMD_WIDTH_HERE == 8
 116 #define CI_TO_CJ_SIMD_4XN(ci)  CI_TO_CJ_J8(ci)
 117 #define X_IND_CI_SIMD_4XN(ci)  X_IND_CI_J8(ci)
 118 #define X_IND_CJ_SIMD_4XN(cj)  X_IND_CJ_J8(cj)
 119 /* Half SIMD with j-cluster size */
 120 #define CI_TO_CJ_SIMD_2XNN(ci) CI_TO_CJ_J4(ci)
 121 #define X_IND_CI_SIMD_2XNN(ci) X_IND_CI_J4(ci)
 122 #define X_IND_CJ_SIMD_2XNN(cj) X_IND_CJ_J4(cj)
 123 #else
 124 #if GMX_SIMD_WIDTH_HERE == 16
 125 #define CI_TO_CJ_SIMD_2XNN(ci) CI_TO_CJ_J8(ci)
 126 #define X_IND_CI_SIMD_2XNN(ci) X_IND_CI_J8(ci)
 127 #define X_IND_CJ_SIMD_2XNN(cj) X_IND_CJ_J8(cj)
 128 #else
 129 #error "unsupported GMX_NBNXN_SIMD_WIDTH"
 130 #endif
 131 #endif
 132 #endif
 133 #endif
 134
 135 #endif /* GMX_NBNXN_SIMD */
 136
 137
 138 #ifdef NBNXN_SEARCH_BB_SIMD4
 139 /* Store bounding boxes corners as quadruplets: xxxxyyyyzzzz */
 140 #define NBNXN_BBXXXX
 141 /* Size of bounding box corners quadruplet */
 142 #define NNBSBB_XXXX      (NNBSBB_D*DIM*STRIDE_PBB)
 143 #endif
 144
 145 /* We shift the i-particles backward for PBC.
 146  * This leads to more conditionals than shifting forward.
 147  * We do this to get more balanced pair lists.
 148  */
 149 #define NBNXN_SHIFT_BACKWARD
 150
 151
 152 /* This define is a lazy way to avoid interdependence of the grid
 153  * and searching data structures.
 154  */
 155 #define NBNXN_NA_SC_MAX (GPU_NSUBCELL*NBNXN_GPU_CLUSTER_SIZE)
 156
 157
 158 static void nbs_cycle_clear(nbnxn_cycle_t *cc)
 159 {
 160     int i;
 161
 162     for (i = 0; i < enbsCCnr; i++)
 163     {
 164         cc[i].count = 0;
 165         cc[i].c     = 0;
 166     }
 167 }
 168
 169 static double Mcyc_av(const nbnxn_cycle_t *cc)
 170 {
 171     return (double)cc->c*1e-6/cc->count;
 172 }
 173
 174 static void nbs_cycle_print(FILE *fp, const nbnxn_search_t nbs)
 175 {
 176     int n;
 177     int t;
 178
 179     fprintf(fp, "\n");
 180     fprintf(fp, "ns %4d grid %4.1f search %4.1f red.f %5.3f",
 181             nbs->cc[enbsCCgrid].count,
 182             Mcyc_av(&nbs->cc[enbsCCgrid]),
 183             Mcyc_av(&nbs->cc[enbsCCsearch]),
 184             Mcyc_av(&nbs->cc[enbsCCreducef]));
 185
 186     if (nbs->nthread_max > 1)
 187     {
 188         if (nbs->cc[enbsCCcombine].count > 0)
 189         {
 190             fprintf(fp, " comb %5.2f",
 191                     Mcyc_av(&nbs->cc[enbsCCcombine]));
 192         }
 193         fprintf(fp, " s. th");
 194         for (t = 0; t < nbs->nthread_max; t++)
 195         {
 196             fprintf(fp, " %4.1f",
 197                     Mcyc_av(&nbs->work[t].cc[enbsCCsearch]));
 198         }
 199     }
 200     fprintf(fp, "\n");
 201 }
 202
 203 static void nbnxn_grid_init(nbnxn_grid_t * grid)
 204 {
 205     grid->cxy_na      = NULL;
 206     grid->cxy_ind     = NULL;
 207     grid->cxy_nalloc  = 0;
 208     grid->bb          = NULL;
 209     grid->bbj         = NULL;
 210     grid->nc_nalloc   = 0;
 211 }
 212
 213 static int get_2log(int n)
 214 {
 215     int log2;
 216
 217     log2 = 0;
 218     while ((1<<log2) < n)
 219     {
 220         log2++;
 221     }
 222     if ((1<<log2) != n)
 223     {
 224         gmx_fatal(FARGS, "nbnxn na_c (%d) is not a power of 2", n);
 225     }
 226
 227     return log2;
 228 }
 229
 230 static int nbnxn_kernel_to_ci_size(int nb_kernel_type)
 231 {
 232     switch (nb_kernel_type)
 233     {
 234         case nbnxnk4x4_PlainC:
 235         case nbnxnk4xN_SIMD_4xN:
 236         case nbnxnk4xN_SIMD_2xNN:
 237             return NBNXN_CPU_CLUSTER_I_SIZE;
 238         case nbnxnk8x8x8_CUDA:
 239         case nbnxnk8x8x8_PlainC:
 240             /* The cluster size for super/sub lists is only set here.
 241              * Any value should work for the pair-search and atomdata code.
 242              * The kernels, of course, might require a particular value.
 243              */
 244             return NBNXN_GPU_CLUSTER_SIZE;
 245         default:
 246             gmx_incons("unknown kernel type");
 247     }
 248
 249     return 0;
 250 }
 251
 252 int nbnxn_kernel_to_cj_size(int nb_kernel_type)
 253 {
 254     int nbnxn_simd_width = 0;
 255     int cj_size          = 0;
 256
 257 #ifdef GMX_NBNXN_SIMD
 258     nbnxn_simd_width = GMX_SIMD_WIDTH_HERE;
 259 #endif
 260
 261     switch (nb_kernel_type)
 262     {
 263         case nbnxnk4x4_PlainC:
 264             cj_size = NBNXN_CPU_CLUSTER_I_SIZE;
 265             break;
 266         case nbnxnk4xN_SIMD_4xN:
 267             cj_size = nbnxn_simd_width;
 268             break;
 269         case nbnxnk4xN_SIMD_2xNN:
 270             cj_size = nbnxn_simd_width/2;
 271             break;
 272         case nbnxnk8x8x8_CUDA:
 273         case nbnxnk8x8x8_PlainC:
 274             cj_size = nbnxn_kernel_to_ci_size(nb_kernel_type);
 275             break;
 276         default:
 277             gmx_incons("unknown kernel type");
 278     }
 279
 280     return cj_size;
 281 }
 282
 283 static int ci_to_cj(int na_cj_2log, int ci)
 284 {
 285     switch (na_cj_2log)
 286     {
 287         case 2: return ci;     break;
 288         case 1: return (ci<<1); break;
 289         case 3: return (ci>>1); break;
 290     }
 291
 292     return 0;
 293 }
 294
 295 gmx_bool nbnxn_kernel_pairlist_simple(int nb_kernel_type)
 296 {
 297     if (nb_kernel_type == nbnxnkNotSet)
 298     {
 299         gmx_fatal(FARGS, "Non-bonded kernel type not set for Verlet-style pair-list.");
 300     }
 301
 302     switch (nb_kernel_type)
 303     {
 304         case nbnxnk8x8x8_CUDA:
 305         case nbnxnk8x8x8_PlainC:
 306             return FALSE;
 307
 308         case nbnxnk4x4_PlainC:
 309         case nbnxnk4xN_SIMD_4xN:
 310         case nbnxnk4xN_SIMD_2xNN:
 311             return TRUE;
 312
 313         default:
 314             gmx_incons("Invalid nonbonded kernel type passed!");
 315             return FALSE;
 316     }
 317 }
 318
 319 void nbnxn_init_search(nbnxn_search_t    * nbs_ptr,
 320                        ivec               *n_dd_cells,
 321                        gmx_domdec_zones_t *zones,
 322                        int                 nthread_max)
 323 {
 324     nbnxn_search_t nbs;
 325     int            d, g, t;
 326
 327     snew(nbs, 1);
 328     *nbs_ptr = nbs;
 329
 330     nbs->DomDec = (n_dd_cells != NULL);
 331
 332     clear_ivec(nbs->dd_dim);
 333     nbs->ngrid = 1;
 334     if (nbs->DomDec)
 335     {
 336         nbs->zones = zones;
 337
 338         for (d = 0; d < DIM; d++)
 339         {
 340             if ((*n_dd_cells)[d] > 1)
 341             {
 342                 nbs->dd_dim[d] = 1;
 343                 /* Each grid matches a DD zone */
 344                 nbs->ngrid *= 2;
 345             }
 346         }
 347     }
 348
 349     snew(nbs->grid, nbs->ngrid);
 350     for (g = 0; g < nbs->ngrid; g++)
 351     {
 352         nbnxn_grid_init(&nbs->grid[g]);
 353     }
 354     nbs->cell        = NULL;
 355     nbs->cell_nalloc = 0;
 356     nbs->a           = NULL;
 357     nbs->a_nalloc    = 0;
 358
 359     nbs->nthread_max = nthread_max;
 360
 361     /* Initialize the work data structures for each thread */
 362     snew(nbs->work, nbs->nthread_max);
 363     for (t = 0; t < nbs->nthread_max; t++)
 364     {
 365         nbs->work[t].cxy_na           = NULL;
 366         nbs->work[t].cxy_na_nalloc    = 0;
 367         nbs->work[t].sort_work        = NULL;
 368         nbs->work[t].sort_work_nalloc = 0;
 369     }
 370
 371     /* Initialize detailed nbsearch cycle counting */
 372     nbs->print_cycles = (getenv("GMX_NBNXN_CYCLE") != 0);
 373     nbs->search_count = 0;
 374     nbs_cycle_clear(nbs->cc);
 375     for (t = 0; t < nbs->nthread_max; t++)
 376     {
 377         nbs_cycle_clear(nbs->work[t].cc);
 378     }
 379 }
 380
 381 static real grid_atom_density(int n, rvec corner0, rvec corner1)
 382 {
 383     rvec size;
 384
 385     rvec_sub(corner1, corner0, size);
 386
 387     return n/(size[XX]*size[YY]*size[ZZ]);
 388 }
 389
 390 static int set_grid_size_xy(const nbnxn_search_t nbs,
 391                             nbnxn_grid_t *grid,
 392                             int dd_zone,
 393                             int n, rvec corner0, rvec corner1,
 394                             real atom_density)
 395 {
 396     rvec size;
 397     int  na_c;
 398     real adens, tlen, tlen_x, tlen_y, nc_max;
 399     int  t;
 400
 401     rvec_sub(corner1, corner0, size);
 402
 403     if (n > grid->na_sc)
 404     {
 405         /* target cell length */
 406         if (grid->bSimple)
 407         {
 408             /* To minimize the zero interactions, we should make
 409              * the largest of the i/j cell cubic.
 410              */
 411             na_c = max(grid->na_c, grid->na_cj);
 412
 413             /* Approximately cubic cells */
 414             tlen   = pow(na_c/atom_density, 1.0/3.0);
 415             tlen_x = tlen;
 416             tlen_y = tlen;
 417         }
 418         else
 419         {
 420             /* Approximately cubic sub cells */
 421             tlen   = pow(grid->na_c/atom_density, 1.0/3.0);
 422             tlen_x = tlen*GPU_NSUBCELL_X;
 423             tlen_y = tlen*GPU_NSUBCELL_Y;
 424         }
 425         /* We round ncx and ncy down, because we get less cell pairs
 426          * in the nbsist when the fixed cell dimensions (x,y) are
 427          * larger than the variable one (z) than the other way around.
 428          */
 429         grid->ncx = max(1, (int)(size[XX]/tlen_x));
 430         grid->ncy = max(1, (int)(size[YY]/tlen_y));
 431     }
 432     else
 433     {
 434         grid->ncx = 1;
 435         grid->ncy = 1;
 436     }
 437
 438     grid->sx     = size[XX]/grid->ncx;
 439     grid->sy     = size[YY]/grid->ncy;
 440     grid->inv_sx = 1/grid->sx;
 441     grid->inv_sy = 1/grid->sy;
 442
 443     if (dd_zone > 0)
 444     {
 445         /* This is a non-home zone, add an extra row of cells
 446          * for particles communicated for bonded interactions.
 447          * These can be beyond the cut-off. It doesn't matter where
 448          * they end up on the grid, but for performance it's better
 449          * if they don't end up in cells that can be within cut-off range.
 450          */
 451         grid->ncx++;
 452         grid->ncy++;
 453     }
 454
 455     /* We need one additional cell entry for particles moved by DD */
 456     if (grid->ncx*grid->ncy+1 > grid->cxy_nalloc)
 457     {
 458         grid->cxy_nalloc = over_alloc_large(grid->ncx*grid->ncy+1);
 459         srenew(grid->cxy_na, grid->cxy_nalloc);
 460         srenew(grid->cxy_ind, grid->cxy_nalloc+1);
 461     }
 462     for (t = 0; t < nbs->nthread_max; t++)
 463     {
 464         if (grid->ncx*grid->ncy+1 > nbs->work[t].cxy_na_nalloc)
 465         {
 466             nbs->work[t].cxy_na_nalloc = over_alloc_large(grid->ncx*grid->ncy+1);
 467             srenew(nbs->work[t].cxy_na, nbs->work[t].cxy_na_nalloc);
 468         }
 469     }
 470
 471     /* Worst case scenario of 1 atom in each last cell */
 472     if (grid->na_cj <= grid->na_c)
 473     {
 474         nc_max = n/grid->na_sc + grid->ncx*grid->ncy;
 475     }
 476     else
 477     {
 478         nc_max = n/grid->na_sc + grid->ncx*grid->ncy*grid->na_cj/grid->na_c;
 479     }
 480
 481     if (nc_max > grid->nc_nalloc)
 482     {
 483         grid->nc_nalloc = over_alloc_large(nc_max);
 484         srenew(grid->nsubc, grid->nc_nalloc);
 485         srenew(grid->bbcz, grid->nc_nalloc*NNBSBB_D);
 486
 487         sfree_aligned(grid->bb);
 488         /* This snew also zeros the contents, this avoid possible
 489          * floating exceptions in SIMD with the unused bb elements.
 490          */
 491         if (grid->bSimple)
 492         {
 493             snew_aligned(grid->bb, grid->nc_nalloc, 16);
 494         }
 495         else
 496         {
 497 #ifdef NBNXN_BBXXXX
 498             int pbb_nalloc;
 499
 500             pbb_nalloc = grid->nc_nalloc*GPU_NSUBCELL/STRIDE_PBB*NNBSBB_XXXX;
 501             snew_aligned(grid->pbb, pbb_nalloc, 16);
 502 #else
 503             snew_aligned(grid->bb, grid->nc_nalloc*GPU_NSUBCELL, 16);
 504 #endif
 505         }
 506
 507         if (grid->bSimple)
 508         {
 509             if (grid->na_cj == grid->na_c)
 510             {
 511                 grid->bbj = grid->bb;
 512             }
 513             else
 514             {
 515                 sfree_aligned(grid->bbj);
 516                 snew_aligned(grid->bbj, grid->nc_nalloc*grid->na_c/grid->na_cj, 16);
 517             }
 518         }
 519
 520         srenew(grid->flags, grid->nc_nalloc);
 521     }
 522
 523     copy_rvec(corner0, grid->c0);
 524     copy_rvec(corner1, grid->c1);
 525
 526     return nc_max;
 527 }
 528
 529 /* We need to sort paricles in grid columns on z-coordinate.
 530  * As particle are very often distributed homogeneously, we a sorting
 531  * algorithm similar to pigeonhole sort. We multiply the z-coordinate
 532  * by a factor, cast to an int and try to store in that hole. If the hole
 533  * is full, we move this or another particle. A second pass is needed to make
 534  * contiguous elements. SORT_GRID_OVERSIZE is the ratio of holes to particles.
 535  * 4 is the optimal value for homogeneous particle distribution and allows
 536  * for an O(#particles) sort up till distributions were all particles are
 537  * concentrated in 1/4 of the space. No NlogN fallback is implemented,
 538  * as it can be expensive to detect imhomogeneous particle distributions.
 539  * SGSF is the maximum ratio of holes used, in the worst case all particles
 540  * end up in the last hole and we need #particles extra holes at the end.
 541  */
 542 #define SORT_GRID_OVERSIZE 4
 543 #define SGSF (SORT_GRID_OVERSIZE + 1)
 544
 545 /* Sort particle index a on coordinates x along dim.
 546  * Backwards tells if we want decreasing iso increasing coordinates.
 547  * h0 is the minimum of the coordinate range.
 548  * invh is the 1/length of the sorting range.
 549  * n_per_h (>=n) is the expected average number of particles per 1/invh
 550  * sort is the sorting work array.
 551  * sort should have a size of at least n_per_h*SORT_GRID_OVERSIZE + n,
 552  * or easier, allocate at least n*SGSF elements.
 553  */
 554 static void sort_atoms(int dim, gmx_bool Backwards,
 555                        int gmx_unused dd_zone,
 556                        int *a, int n, rvec *x,
 557                        real h0, real invh, int n_per_h,
 558                        int *sort)
 559 {
 560     int nsort, i, c;
 561     int zi, zim, zi_min, zi_max;
 562     int cp, tmp;
 563
 564     if (n <= 1)
 565     {
 566         /* Nothing to do */
 567         return;
 568     }
 569
 570 #ifndef NDEBUG
 571     if (n > n_per_h)
 572     {
 573         gmx_incons("n > n_per_h");
 574     }
 575 #endif
 576
 577     /* Transform the inverse range height into the inverse hole height */
 578     invh *= n_per_h*SORT_GRID_OVERSIZE;
 579
 580     /* Set nsort to the maximum possible number of holes used.
 581      * In worst case all n elements end up in the last bin.
 582      */
 583     nsort = n_per_h*SORT_GRID_OVERSIZE + n;
 584
 585     /* Determine the index range used, so we can limit it for the second pass */
 586     zi_min = INT_MAX;
 587     zi_max = -1;
 588
 589     /* Sort the particles using a simple index sort */
 590     for (i = 0; i < n; i++)
 591     {
 592         /* The cast takes care of float-point rounding effects below zero.
 593          * This code assumes particles are less than 1/SORT_GRID_OVERSIZE
 594          * times the box height out of the box.
 595          */
 596         zi = (int)((x[a[i]][dim] - h0)*invh);
 597
 598 #ifndef NDEBUG
 599         /* As we can have rounding effect, we use > iso >= here */
 600         if (zi < 0 || (dd_zone == 0 && zi > n_per_h*SORT_GRID_OVERSIZE))
 601         {
 602             gmx_fatal(FARGS, "(int)((x[%d][%c]=%f - %f)*%f) = %d, not in 0 - %d*%d\n",
 603                       a[i], 'x'+dim, x[a[i]][dim], h0, invh, zi,
 604                       n_per_h, SORT_GRID_OVERSIZE);
 605         }
 606 #endif
 607
 608         /* In a non-local domain, particles communcated for bonded interactions
 609          * can be far beyond the grid size, which is set by the non-bonded
 610          * cut-off distance. We sort such particles into the last cell.
 611          */
 612         if (zi > n_per_h*SORT_GRID_OVERSIZE)
 613         {
 614             zi = n_per_h*SORT_GRID_OVERSIZE;
 615         }
 616
 617         /* Ideally this particle should go in sort cell zi,
 618          * but that might already be in use,
 619          * in that case find the first empty cell higher up
 620          */
 621         if (sort[zi] < 0)
 622         {
 623             sort[zi] = a[i];
 624             zi_min   = min(zi_min, zi);
 625             zi_max   = max(zi_max, zi);
 626         }
 627         else
 628         {
 629             /* We have multiple atoms in the same sorting slot.
 630              * Sort on real z for minimal bounding box size.
 631              * There is an extra check for identical z to ensure
 632              * well-defined output order, independent of input order
 633              * to ensure binary reproducibility after restarts.
 634              */
 635             while (sort[zi] >= 0 && ( x[a[i]][dim] >  x[sort[zi]][dim] ||
 636                                       (x[a[i]][dim] == x[sort[zi]][dim] &&
 637                                        a[i] > sort[zi])))
 638             {
 639                 zi++;
 640             }
 641
 642             if (sort[zi] >= 0)
 643             {
 644                 /* Shift all elements by one slot until we find an empty slot */
 645                 cp  = sort[zi];
 646                 zim = zi + 1;
 647                 while (sort[zim] >= 0)
 648                 {
 649                     tmp       = sort[zim];
 650                     sort[zim] = cp;
 651                     cp        = tmp;
 652                     zim++;
 653                 }
 654                 sort[zim] = cp;
 655                 zi_max    = max(zi_max, zim);
 656             }
 657             sort[zi] = a[i];
 658             zi_max   = max(zi_max, zi);
 659         }
 660     }
 661
 662     c = 0;
 663     if (!Backwards)
 664     {
 665         for (zi = 0; zi < nsort; zi++)
 666         {
 667             if (sort[zi] >= 0)
 668             {
 669                 a[c++]   = sort[zi];
 670                 sort[zi] = -1;
 671             }
 672         }
 673     }
 674     else
 675     {
 676         for (zi = zi_max; zi >= zi_min; zi--)
 677         {
 678             if (sort[zi] >= 0)
 679             {
 680                 a[c++]   = sort[zi];
 681                 sort[zi] = -1;
 682             }
 683         }
 684     }
 685     if (c < n)
 686     {
 687         gmx_incons("Lost particles while sorting");
 688     }
 689 }
 690
 691 #ifdef GMX_DOUBLE
 692 #define R2F_D(x) ((float)((x) >= 0 ? ((1-GMX_FLOAT_EPS)*(x)) : ((1+GMX_FLOAT_EPS)*(x))))
 693 #define R2F_U(x) ((float)((x) >= 0 ? ((1+GMX_FLOAT_EPS)*(x)) : ((1-GMX_FLOAT_EPS)*(x))))
 694 #else
 695 #define R2F_D(x) (x)
 696 #define R2F_U(x) (x)
 697 #endif
 698
 699 /* Coordinate order x,y,z, bb order xyz0 */
 700 static void calc_bounding_box(int na, int stride, const real *x, nbnxn_bb_t *bb)
 701 {
 702     int  i, j;
 703     real xl, xh, yl, yh, zl, zh;
 704
 705     i  = 0;
 706     xl = x[i+XX];
 707     xh = x[i+XX];
 708     yl = x[i+YY];
 709     yh = x[i+YY];
 710     zl = x[i+ZZ];
 711     zh = x[i+ZZ];
 712     i += stride;
 713     for (j = 1; j < na; j++)
 714     {
 715         xl = min(xl, x[i+XX]);
 716         xh = max(xh, x[i+XX]);
 717         yl = min(yl, x[i+YY]);
 718         yh = max(yh, x[i+YY]);
 719         zl = min(zl, x[i+ZZ]);
 720         zh = max(zh, x[i+ZZ]);
 721         i += stride;
 722     }
 723     /* Note: possible double to float conversion here */
 724     bb->lower[BB_X] = R2F_D(xl);
 725     bb->lower[BB_Y] = R2F_D(yl);
 726     bb->lower[BB_Z] = R2F_D(zl);
 727     bb->upper[BB_X] = R2F_U(xh);
 728     bb->upper[BB_Y] = R2F_U(yh);
 729     bb->upper[BB_Z] = R2F_U(zh);
 730 }
 731
 732 /* Packed coordinates, bb order xyz0 */
 733 static void calc_bounding_box_x_x4(int na, const real *x, nbnxn_bb_t *bb)
 734 {
 735     int  j;
 736     real xl, xh, yl, yh, zl, zh;
 737
 738     xl = x[XX*PACK_X4];
 739     xh = x[XX*PACK_X4];
 740     yl = x[YY*PACK_X4];
 741     yh = x[YY*PACK_X4];
 742     zl = x[ZZ*PACK_X4];
 743     zh = x[ZZ*PACK_X4];
 744     for (j = 1; j < na; j++)
 745     {
 746         xl = min(xl, x[j+XX*PACK_X4]);
 747         xh = max(xh, x[j+XX*PACK_X4]);
 748         yl = min(yl, x[j+YY*PACK_X4]);
 749         yh = max(yh, x[j+YY*PACK_X4]);
 750         zl = min(zl, x[j+ZZ*PACK_X4]);
 751         zh = max(zh, x[j+ZZ*PACK_X4]);
 752     }
 753     /* Note: possible double to float conversion here */
 754     bb->lower[BB_X] = R2F_D(xl);
 755     bb->lower[BB_Y] = R2F_D(yl);
 756     bb->lower[BB_Z] = R2F_D(zl);
 757     bb->upper[BB_X] = R2F_U(xh);
 758     bb->upper[BB_Y] = R2F_U(yh);
 759     bb->upper[BB_Z] = R2F_U(zh);
 760 }
 761
 762 /* Packed coordinates, bb order xyz0 */
 763 static void calc_bounding_box_x_x8(int na, const real *x, nbnxn_bb_t *bb)
 764 {
 765     int  j;
 766     real xl, xh, yl, yh, zl, zh;
 767
 768     xl = x[XX*PACK_X8];
 769     xh = x[XX*PACK_X8];
 770     yl = x[YY*PACK_X8];
 771     yh = x[YY*PACK_X8];
 772     zl = x[ZZ*PACK_X8];
 773     zh = x[ZZ*PACK_X8];
 774     for (j = 1; j < na; j++)
 775     {
 776         xl = min(xl, x[j+XX*PACK_X8]);
 777         xh = max(xh, x[j+XX*PACK_X8]);
 778         yl = min(yl, x[j+YY*PACK_X8]);
 779         yh = max(yh, x[j+YY*PACK_X8]);
 780         zl = min(zl, x[j+ZZ*PACK_X8]);
 781         zh = max(zh, x[j+ZZ*PACK_X8]);
 782     }
 783     /* Note: possible double to float conversion here */
 784     bb->lower[BB_X] = R2F_D(xl);
 785     bb->lower[BB_Y] = R2F_D(yl);
 786     bb->lower[BB_Z] = R2F_D(zl);
 787     bb->upper[BB_X] = R2F_U(xh);
 788     bb->upper[BB_Y] = R2F_U(yh);
 789     bb->upper[BB_Z] = R2F_U(zh);
 790 }
 791
 792 /* Packed coordinates, bb order xyz0 */
 793 static void calc_bounding_box_x_x4_halves(int na, const real *x,
 794                                           nbnxn_bb_t *bb, nbnxn_bb_t *bbj)
 795 {
 796     calc_bounding_box_x_x4(min(na, 2), x, bbj);
 797
 798     if (na > 2)
 799     {
 800         calc_bounding_box_x_x4(min(na-2, 2), x+(PACK_X4>>1), bbj+1);
 801     }
 802     else
 803     {
 804         /* Set the "empty" bounding box to the same as the first one,
 805          * so we don't need to treat special cases in the rest of the code.
 806          */
 807 #ifdef NBNXN_SEARCH_BB_SIMD4
 808         gmx_simd4_store_pr(&bbj[1].lower[0], gmx_simd4_load_bb_pr(&bbj[0].lower[0]));
 809         gmx_simd4_store_pr(&bbj[1].upper[0], gmx_simd4_load_bb_pr(&bbj[0].upper[0]));
 810 #else
 811         bbj[1] = bbj[0];
 812 #endif
 813     }
 814
 815 #ifdef NBNXN_SEARCH_BB_SIMD4
 816     gmx_simd4_store_pr(&bb->lower[0],
 817                        gmx_simd4_min_pr(gmx_simd4_load_bb_pr(&bbj[0].lower[0]),
 818                                         gmx_simd4_load_bb_pr(&bbj[1].lower[0])));
 819     gmx_simd4_store_pr(&bb->upper[0],
 820                        gmx_simd4_max_pr(gmx_simd4_load_bb_pr(&bbj[0].upper[0]),
 821                                         gmx_simd4_load_bb_pr(&bbj[1].upper[0])));
 822 #else
 823     {
 824         int i;
 825
 826         for (i = 0; i < NNBSBB_C; i++)
 827         {
 828             bb->lower[i] = min(bbj[0].lower[i], bbj[1].lower[i]);
 829             bb->upper[i] = max(bbj[0].upper[i], bbj[1].upper[i]);
 830         }
 831     }
 832 #endif
 833 }
 834
 835 #ifdef NBNXN_SEARCH_BB_SIMD4
 836
 837 /* Coordinate order xyz, bb order xxxxyyyyzzzz */
 838 static void calc_bounding_box_xxxx(int na, int stride, const real *x, float *bb)
 839 {
 840     int  i, j;
 841     real xl, xh, yl, yh, zl, zh;
 842
 843     i  = 0;
 844     xl = x[i+XX];
 845     xh = x[i+XX];
 846     yl = x[i+YY];
 847     yh = x[i+YY];
 848     zl = x[i+ZZ];
 849     zh = x[i+ZZ];
 850     i += stride;
 851     for (j = 1; j < na; j++)
 852     {
 853         xl = min(xl, x[i+XX]);
 854         xh = max(xh, x[i+XX]);
 855         yl = min(yl, x[i+YY]);
 856         yh = max(yh, x[i+YY]);
 857         zl = min(zl, x[i+ZZ]);
 858         zh = max(zh, x[i+ZZ]);
 859         i += stride;
 860     }
 861     /* Note: possible double to float conversion here */
 862     bb[0*STRIDE_PBB] = R2F_D(xl);
 863     bb[1*STRIDE_PBB] = R2F_D(yl);
 864     bb[2*STRIDE_PBB] = R2F_D(zl);
 865     bb[3*STRIDE_PBB] = R2F_U(xh);
 866     bb[4*STRIDE_PBB] = R2F_U(yh);
 867     bb[5*STRIDE_PBB] = R2F_U(zh);
 868 }
 869
 870 #endif /* NBNXN_SEARCH_BB_SIMD4 */
 871
 872 #ifdef NBNXN_SEARCH_SIMD4_FLOAT_X_BB
 873
 874 /* Coordinate order xyz?, bb order xyz0 */
 875 static void calc_bounding_box_simd4(int na, const float *x, nbnxn_bb_t *bb)
 876 {
 877     gmx_simd4_pr bb_0_S, bb_1_S;
 878     gmx_simd4_pr x_S;
 879
 880     int          i;
 881
 882     bb_0_S = gmx_simd4_load_bb_pr(x);
 883     bb_1_S = bb_0_S;
 884
 885     for (i = 1; i < na; i++)
 886     {
 887         x_S    = gmx_simd4_load_bb_pr(x+i*NNBSBB_C);
 888         bb_0_S = gmx_simd4_min_pr(bb_0_S, x_S);
 889         bb_1_S = gmx_simd4_max_pr(bb_1_S, x_S);
 890     }
 891
 892     gmx_simd4_store_pr(&bb->lower[0], bb_0_S);
 893     gmx_simd4_store_pr(&bb->upper[0], bb_1_S);
 894 }
 895
 896 /* Coordinate order xyz?, bb order xxxxyyyyzzzz */
 897 static void calc_bounding_box_xxxx_simd4(int na, const float *x,
 898                                          nbnxn_bb_t *bb_work_aligned,
 899                                          real *bb)
 900 {
 901     calc_bounding_box_simd4(na, x, bb_work_aligned);
 902
 903     bb[0*STRIDE_PBB] = bb_work_aligned->lower[BB_X];
 904     bb[1*STRIDE_PBB] = bb_work_aligned->lower[BB_Y];
 905     bb[2*STRIDE_PBB] = bb_work_aligned->lower[BB_Z];
 906     bb[3*STRIDE_PBB] = bb_work_aligned->upper[BB_X];
 907     bb[4*STRIDE_PBB] = bb_work_aligned->upper[BB_Y];
 908     bb[5*STRIDE_PBB] = bb_work_aligned->upper[BB_Z];
 909 }
 910
 911 #endif /* NBNXN_SEARCH_SIMD4_FLOAT_X_BB */
 912
 913
 914 /* Combines pairs of consecutive bounding boxes */
 915 static void combine_bounding_box_pairs(nbnxn_grid_t *grid, const nbnxn_bb_t *bb)
 916 {
 917     int    i, j, sc2, nc2, c2;
 918
 919     for (i = 0; i < grid->ncx*grid->ncy; i++)
 920     {
 921         /* Starting bb in a column is expected to be 2-aligned */
 922         sc2 = grid->cxy_ind[i]>>1;
 923         /* For odd numbers skip the last bb here */
 924         nc2 = (grid->cxy_na[i]+3)>>(2+1);
 925         for (c2 = sc2; c2 < sc2+nc2; c2++)
 926         {
 927 #ifdef NBNXN_SEARCH_BB_SIMD4
 928             gmx_simd4_pr min_S, max_S;
 929
 930             min_S = gmx_simd4_min_pr(gmx_simd4_load_bb_pr(&bb[c2*2+0].lower[0]),
 931                                      gmx_simd4_load_bb_pr(&bb[c2*2+1].lower[0]));
 932             max_S = gmx_simd4_max_pr(gmx_simd4_load_bb_pr(&bb[c2*2+0].upper[0]),
 933                                      gmx_simd4_load_bb_pr(&bb[c2*2+1].upper[0]));
 934             gmx_simd4_store_pr(&grid->bbj[c2].lower[0], min_S);
 935             gmx_simd4_store_pr(&grid->bbj[c2].upper[0], max_S);
 936 #else
 937             for (j = 0; j < NNBSBB_C; j++)
 938             {
 939                 grid->bbj[c2].lower[j] = min(bb[c2*2+0].lower[j],
 940                                              bb[c2*2+1].lower[j]);
 941                 grid->bbj[c2].upper[j] = max(bb[c2*2+0].upper[j],
 942                                              bb[c2*2+1].upper[j]);
 943             }
 944 #endif
 945         }
 946         if (((grid->cxy_na[i]+3)>>2) & 1)
 947         {
 948             /* The bb count in this column is odd: duplicate the last bb */
 949             for (j = 0; j < NNBSBB_C; j++)
 950             {
 951                 grid->bbj[c2].lower[j] = bb[c2*2].lower[j];
 952                 grid->bbj[c2].upper[j] = bb[c2*2].upper[j];
 953             }
 954         }
 955     }
 956 }
 957
 958
 959 /* Prints the average bb size, used for debug output */
 960 static void print_bbsizes_simple(FILE                *fp,
 961                                  const nbnxn_search_t nbs,
 962                                  const nbnxn_grid_t  *grid)
 963 {
 964     int  c, d;
 965     dvec ba;
 966
 967     clear_dvec(ba);
 968     for (c = 0; c < grid->nc; c++)
 969     {
 970         for (d = 0; d < DIM; d++)
 971         {
 972             ba[d] += grid->bb[c].upper[d] - grid->bb[c].lower[d];
 973         }
 974     }
 975     dsvmul(1.0/grid->nc, ba, ba);
 976
 977     fprintf(fp, "ns bb: %4.2f %4.2f %4.2f  %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
 978             nbs->box[XX][XX]/grid->ncx,
 979             nbs->box[YY][YY]/grid->ncy,
 980             nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/grid->nc,
 981             ba[XX], ba[YY], ba[ZZ],
 982             ba[XX]*grid->ncx/nbs->box[XX][XX],
 983             ba[YY]*grid->ncy/nbs->box[YY][YY],
 984             ba[ZZ]*grid->nc/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
 985 }
 986
 987 /* Prints the average bb size, used for debug output */
 988 static void print_bbsizes_supersub(FILE                *fp,
 989                                    const nbnxn_search_t nbs,
 990                                    const nbnxn_grid_t  *grid)
 991 {
 992     int  ns, c, s;
 993     dvec ba;
 994
 995     clear_dvec(ba);
 996     ns = 0;
 997     for (c = 0; c < grid->nc; c++)
 998     {
 999 #ifdef NBNXN_BBXXXX
1000         for (s = 0; s < grid->nsubc[c]; s += STRIDE_PBB)
1001         {
1002             int cs_w, i, d;
1003
1004             cs_w = (c*GPU_NSUBCELL + s)/STRIDE_PBB;
1005             for (i = 0; i < STRIDE_PBB; i++)
1006             {
1007                 for (d = 0; d < DIM; d++)
1008                 {
1009                     ba[d] +=
1010                         grid->pbb[cs_w*NNBSBB_XXXX+(DIM+d)*STRIDE_PBB+i] -
1011                         grid->pbb[cs_w*NNBSBB_XXXX+     d *STRIDE_PBB+i];
1012                 }
1013             }
1014         }
1015 #else
1016         for (s = 0; s < grid->nsubc[c]; s++)
1017         {
1018             int cs, d;
1019
1020             cs = c*GPU_NSUBCELL + s;
1021             for (d = 0; d < DIM; d++)
1022             {
1023                 ba[d] += grid->bb[cs].upper[d] - grid->bb[cs].lower[d];
1024             }
1025         }
1026 #endif
1027         ns += grid->nsubc[c];
1028     }
1029     dsvmul(1.0/ns, ba, ba);
1030
1031     fprintf(fp, "ns bb: %4.2f %4.2f %4.2f  %4.2f %4.2f %4.2f rel %4.2f %4.2f %4.2f\n",
1032             nbs->box[XX][XX]/(grid->ncx*GPU_NSUBCELL_X),
1033             nbs->box[YY][YY]/(grid->ncy*GPU_NSUBCELL_Y),
1034             nbs->box[ZZ][ZZ]*grid->ncx*grid->ncy/(grid->nc*GPU_NSUBCELL_Z),
1035             ba[XX], ba[YY], ba[ZZ],
1036             ba[XX]*grid->ncx*GPU_NSUBCELL_X/nbs->box[XX][XX],
1037             ba[YY]*grid->ncy*GPU_NSUBCELL_Y/nbs->box[YY][YY],
1038             ba[ZZ]*grid->nc*GPU_NSUBCELL_Z/(grid->ncx*grid->ncy*nbs->box[ZZ][ZZ]));
1039 }
1040
1041 /* Potentially sorts atoms on LJ coefficients !=0 and ==0.
1042  * Also sets interaction flags.
1043  */
1044 void sort_on_lj(int na_c,
1045                 int a0, int a1, const int *atinfo,
1046                 int *order,
1047                 int *flags)
1048 {
1049     int      subc, s, a, n1, n2, a_lj_max, i, j;
1050     int      sort1[NBNXN_NA_SC_MAX/GPU_NSUBCELL];
1051     int      sort2[NBNXN_NA_SC_MAX/GPU_NSUBCELL];
1052     gmx_bool haveQ;
1053
1054     *flags = 0;
1055
1056     subc = 0;
1057     for (s = a0; s < a1; s += na_c)
1058     {
1059         /* Make lists for this (sub-)cell on atoms with and without LJ */
1060         n1       = 0;
1061         n2       = 0;
1062         haveQ    = FALSE;
1063         a_lj_max = -1;
1064         for (a = s; a < min(s+na_c, a1); a++)
1065         {
1066             haveQ = haveQ || GET_CGINFO_HAS_Q(atinfo[order[a]]);
1067
1068             if (GET_CGINFO_HAS_VDW(atinfo[order[a]]))
1069             {
1070                 sort1[n1++] = order[a];
1071                 a_lj_max    = a;
1072             }
1073             else
1074             {
1075                 sort2[n2++] = order[a];
1076             }
1077         }
1078
1079         /* If we don't have atom with LJ, there's nothing to sort */
1080         if (n1 > 0)
1081         {
1082             *flags |= NBNXN_CI_DO_LJ(subc);
1083
1084             if (2*n1 <= na_c)
1085             {
1086                 /* Only sort when strictly necessary. Ordering particles
1087                  * Ordering particles can lead to less accurate summation
1088                  * due to rounding, both for LJ and Coulomb interactions.
1089                  */
1090                 if (2*(a_lj_max - s) >= na_c)
1091                 {
1092                     for (i = 0; i < n1; i++)
1093                     {
1094                         order[a0+i] = sort1[i];
1095                     }
1096                     for (j = 0; j < n2; j++)
1097                     {
1098                         order[a0+n1+j] = sort2[j];
1099                     }
1100                 }
1101
1102                 *flags |= NBNXN_CI_HALF_LJ(subc);
1103             }
1104         }
1105         if (haveQ)
1106         {
1107             *flags |= NBNXN_CI_DO_COUL(subc);
1108         }
1109         subc++;
1110     }
1111 }
1112
1113 /* Fill a pair search cell with atoms.
1114  * Potentially sorts atoms and sets the interaction flags.
1115  */
1116 void fill_cell(const nbnxn_search_t nbs,
1117                nbnxn_grid_t *grid,
1118                nbnxn_atomdata_t *nbat,
1119                int a0, int a1,
1120                const int *atinfo,
1121                rvec *x,
1122                int sx, int sy, int sz,
1123                nbnxn_bb_t gmx_unused *bb_work_aligned)
1124 {
1125     int         na, a;
1126     size_t      offset;
1127     nbnxn_bb_t *bb_ptr;
1128 #ifdef NBNXN_BBXXXX
1129     float      *pbb_ptr;
1130 #endif
1131
1132     na = a1 - a0;
1133
1134     if (grid->bSimple)
1135     {
1136         sort_on_lj(grid->na_c, a0, a1, atinfo, nbs->a,
1137                    grid->flags+(a0>>grid->na_c_2log)-grid->cell0);
1138     }
1139
1140     /* Now we have sorted the atoms, set the cell indices */
1141     for (a = a0; a < a1; a++)
1142     {
1143         nbs->cell[nbs->a[a]] = a;
1144     }
1145
1146     copy_rvec_to_nbat_real(nbs->a+a0, a1-a0, grid->na_c, x,
1147                            nbat->XFormat, nbat->x, a0,
1148                            sx, sy, sz);
1149
1150     if (nbat->XFormat == nbatX4)
1151     {
1152         /* Store the bounding boxes as xyz.xyz. */
1153         offset = (a0 - grid->cell0*grid->na_sc) >> grid->na_c_2log;
1154         bb_ptr = grid->bb + offset;
1155
1156 #if defined GMX_NBNXN_SIMD && GMX_SIMD_WIDTH_HERE == 2
1157         if (2*grid->na_cj == grid->na_c)
1158         {
1159             calc_bounding_box_x_x4_halves(na, nbat->x+X4_IND_A(a0), bb_ptr,
1160                                           grid->bbj+offset*2);
1161         }
1162         else
1163 #endif
1164         {
1165             calc_bounding_box_x_x4(na, nbat->x+X4_IND_A(a0), bb_ptr);
1166         }
1167     }
1168     else if (nbat->XFormat == nbatX8)
1169     {
1170         /* Store the bounding boxes as xyz.xyz. */
1171         offset = (a0 - grid->cell0*grid->na_sc) >> grid->na_c_2log;
1172         bb_ptr = grid->bb + offset;
1173
1174         calc_bounding_box_x_x8(na, nbat->x+X8_IND_A(a0), bb_ptr);
1175     }
1176 #ifdef NBNXN_BBXXXX
1177     else if (!grid->bSimple)
1178     {
1179         /* Store the bounding boxes in a format convenient
1180          * for SIMD4 calculations: xxxxyyyyzzzz...
1181          */
1182         pbb_ptr =
1183             grid->pbb +
1184             ((a0-grid->cell0*grid->na_sc)>>(grid->na_c_2log+STRIDE_PBB_2LOG))*NNBSBB_XXXX +
1185             (((a0-grid->cell0*grid->na_sc)>>grid->na_c_2log) & (STRIDE_PBB-1));
1186
1187 #ifdef NBNXN_SEARCH_SIMD4_FLOAT_X_BB
1188         if (nbat->XFormat == nbatXYZQ)
1189         {
1190             calc_bounding_box_xxxx_simd4(na, nbat->x+a0*nbat->xstride,
1191                                          bb_work_aligned, pbb_ptr);
1192         }
1193         else
1194 #endif
1195         {
1196             calc_bounding_box_xxxx(na, nbat->xstride, nbat->x+a0*nbat->xstride,
1197                                    pbb_ptr);
1198         }
1199         if (gmx_debug_at)
1200         {
1201             fprintf(debug, "%2d %2d %2d bb %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
1202                     sx, sy, sz,
1203                     pbb_ptr[0*STRIDE_PBB], pbb_ptr[3*STRIDE_PBB],
1204                     pbb_ptr[1*STRIDE_PBB], pbb_ptr[4*STRIDE_PBB],
1205                     pbb_ptr[2*STRIDE_PBB], pbb_ptr[5*STRIDE_PBB]);
1206         }
1207     }
1208 #endif
1209     else
1210     {
1211         /* Store the bounding boxes as xyz.xyz. */
1212         bb_ptr = grid->bb+((a0-grid->cell0*grid->na_sc)>>grid->na_c_2log);
1213
1214         calc_bounding_box(na, nbat->xstride, nbat->x+a0*nbat->xstride,
1215                           bb_ptr);
1216
1217         if (gmx_debug_at)
1218         {
1219             int bbo;
1220             bbo = (a0 - grid->cell0*grid->na_sc)/grid->na_c;
1221             fprintf(debug, "%2d %2d %2d bb %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f\n",
1222                     sx, sy, sz,
1223                     grid->bb[bbo].lower[BB_X],
1224                     grid->bb[bbo].lower[BB_Y],
1225                     grid->bb[bbo].lower[BB_Z],
1226                     grid->bb[bbo].upper[BB_X],
1227                     grid->bb[bbo].upper[BB_Y],
1228                     grid->bb[bbo].upper[BB_Z]);
1229         }
1230     }
1231 }
1232
1233 /* Spatially sort the atoms within one grid column */
1234 static void sort_columns_simple(const nbnxn_search_t nbs,
1235                                 int dd_zone,
1236                                 nbnxn_grid_t *grid,
1237                                 int a0, int a1,
1238                                 const int *atinfo,
1239                                 rvec *x,
1240                                 nbnxn_atomdata_t *nbat,
1241                                 int cxy_start, int cxy_end,
1242                                 int *sort_work)
1243 {
1244     int  cxy;
1245     int  cx, cy, cz, ncz, cfilled, c;
1246     int  na, ash, ind, a;
1247     int  na_c, ash_c;
1248
1249     if (debug)
1250     {
1251         fprintf(debug, "cell0 %d sorting columns %d - %d, atoms %d - %d\n",
1252                 grid->cell0, cxy_start, cxy_end, a0, a1);
1253     }
1254
1255     /* Sort the atoms within each x,y column in 3 dimensions */
1256     for (cxy = cxy_start; cxy < cxy_end; cxy++)
1257     {
1258         cx = cxy/grid->ncy;
1259         cy = cxy - cx*grid->ncy;
1260
1261         na  = grid->cxy_na[cxy];
1262         ncz = grid->cxy_ind[cxy+1] - grid->cxy_ind[cxy];
1263         ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
1264
1265         /* Sort the atoms within each x,y column on z coordinate */
1266         sort_atoms(ZZ, FALSE, dd_zone,
1267                    nbs->a+ash, na, x,
1268                    grid->c0[ZZ],
1269                    1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
1270                    sort_work);
1271
1272         /* Fill the ncz cells in this column */
1273         cfilled = grid->cxy_ind[cxy];
1274         for (cz = 0; cz < ncz; cz++)
1275         {
1276             c  = grid->cxy_ind[cxy] + cz;
1277
1278             ash_c = ash + cz*grid->na_sc;
1279             na_c  = min(grid->na_sc, na-(ash_c-ash));
1280
1281             fill_cell(nbs, grid, nbat,
1282                       ash_c, ash_c+na_c, atinfo, x,
1283                       grid->na_sc*cx + (dd_zone >> 2),
1284                       grid->na_sc*cy + (dd_zone & 3),
1285                       grid->na_sc*cz,
1286                       NULL);
1287
1288             /* This copy to bbcz is not really necessary.
1289              * But it allows to use the same grid search code
1290              * for the simple and supersub cell setups.
1291              */
1292             if (na_c > 0)
1293             {
1294                 cfilled = c;
1295             }
1296             grid->bbcz[c*NNBSBB_D  ] = grid->bb[cfilled].lower[BB_Z];
1297             grid->bbcz[c*NNBSBB_D+1] = grid->bb[cfilled].upper[BB_Z];
1298         }
1299
1300         /* Set the unused atom indices to -1 */
1301         for (ind = na; ind < ncz*grid->na_sc; ind++)
1302         {
1303             nbs->a[ash+ind] = -1;
1304         }
1305     }
1306 }
1307
1308 /* Spatially sort the atoms within one grid column */
1309 static void sort_columns_supersub(const nbnxn_search_t nbs,
1310                                   int dd_zone,
1311                                   nbnxn_grid_t *grid,
1312                                   int a0, int a1,
1313                                   const int *atinfo,
1314                                   rvec *x,
1315                                   nbnxn_atomdata_t *nbat,
1316                                   int cxy_start, int cxy_end,
1317                                   int *sort_work)
1318 {
1319     int  cxy;
1320     int  cx, cy, cz = -1, c = -1, ncz;
1321     int  na, ash, na_c, ind, a;
1322     int  subdiv_z, sub_z, na_z, ash_z;
1323     int  subdiv_y, sub_y, na_y, ash_y;
1324     int  subdiv_x, sub_x, na_x, ash_x;
1325
1326     /* cppcheck-suppress unassignedVariable */
1327     nbnxn_bb_t bb_work_array[2], *bb_work_aligned;
1328
1329     bb_work_aligned = (nbnxn_bb_t *)(((size_t)(bb_work_array+1)) & (~((size_t)15)));
1330
1331     if (debug)
1332     {
1333         fprintf(debug, "cell0 %d sorting columns %d - %d, atoms %d - %d\n",
1334                 grid->cell0, cxy_start, cxy_end, a0, a1);
1335     }
1336
1337     subdiv_x = grid->na_c;
1338     subdiv_y = GPU_NSUBCELL_X*subdiv_x;
1339     subdiv_z = GPU_NSUBCELL_Y*subdiv_y;
1340
1341     /* Sort the atoms within each x,y column in 3 dimensions */
1342     for (cxy = cxy_start; cxy < cxy_end; cxy++)
1343     {
1344         cx = cxy/grid->ncy;
1345         cy = cxy - cx*grid->ncy;
1346
1347         na  = grid->cxy_na[cxy];
1348         ncz = grid->cxy_ind[cxy+1] - grid->cxy_ind[cxy];
1349         ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
1350
1351         /* Sort the atoms within each x,y column on z coordinate */
1352         sort_atoms(ZZ, FALSE, dd_zone,
1353                    nbs->a+ash, na, x,
1354                    grid->c0[ZZ],
1355                    1.0/nbs->box[ZZ][ZZ], ncz*grid->na_sc,
1356                    sort_work);
1357
1358         /* This loop goes over the supercells and subcells along z at once */
1359         for (sub_z = 0; sub_z < ncz*GPU_NSUBCELL_Z; sub_z++)
1360         {
1361             ash_z = ash + sub_z*subdiv_z;
1362             na_z  = min(subdiv_z, na-(ash_z-ash));
1363
1364             /* We have already sorted on z */
1365
1366             if (sub_z % GPU_NSUBCELL_Z == 0)
1367             {
1368                 cz = sub_z/GPU_NSUBCELL_Z;
1369                 c  = grid->cxy_ind[cxy] + cz;
1370
1371                 /* The number of atoms in this supercell */
1372                 na_c = min(grid->na_sc, na-(ash_z-ash));
1373
1374                 grid->nsubc[c] = min(GPU_NSUBCELL, (na_c+grid->na_c-1)/grid->na_c);
1375
1376                 /* Store the z-boundaries of the super cell */
1377                 grid->bbcz[c*NNBSBB_D  ] = x[nbs->a[ash_z]][ZZ];
1378                 grid->bbcz[c*NNBSBB_D+1] = x[nbs->a[ash_z+na_c-1]][ZZ];
1379             }
1380
1381 #if GPU_NSUBCELL_Y > 1
1382             /* Sort the atoms along y */
1383             sort_atoms(YY, (sub_z & 1), dd_zone,
1384                        nbs->a+ash_z, na_z, x,
1385                        grid->c0[YY]+cy*grid->sy,
1386                        grid->inv_sy, subdiv_z,
1387                        sort_work);
1388 #endif
1389
1390             for (sub_y = 0; sub_y < GPU_NSUBCELL_Y; sub_y++)
1391             {
1392                 ash_y = ash_z + sub_y*subdiv_y;
1393                 na_y  = min(subdiv_y, na-(ash_y-ash));
1394
1395 #if GPU_NSUBCELL_X > 1
1396                 /* Sort the atoms along x */
1397                 sort_atoms(XX, ((cz*GPU_NSUBCELL_Y + sub_y) & 1), dd_zone,
1398                            nbs->a+ash_y, na_y, x,
1399                            grid->c0[XX]+cx*grid->sx,
1400                            grid->inv_sx, subdiv_y,
1401                            sort_work);
1402 #endif
1403
1404                 for (sub_x = 0; sub_x < GPU_NSUBCELL_X; sub_x++)
1405                 {
1406                     ash_x = ash_y + sub_x*subdiv_x;
1407                     na_x  = min(subdiv_x, na-(ash_x-ash));
1408
1409                     fill_cell(nbs, grid, nbat,
1410                               ash_x, ash_x+na_x, atinfo, x,
1411                               grid->na_c*(cx*GPU_NSUBCELL_X+sub_x) + (dd_zone >> 2),
1412                               grid->na_c*(cy*GPU_NSUBCELL_Y+sub_y) + (dd_zone & 3),
1413                               grid->na_c*sub_z,
1414                               bb_work_aligned);
1415                 }
1416             }
1417         }
1418
1419         /* Set the unused atom indices to -1 */
1420         for (ind = na; ind < ncz*grid->na_sc; ind++)
1421         {
1422             nbs->a[ash+ind] = -1;
1423         }
1424     }
1425 }
1426
1427 /* Determine in which grid column atoms should go */
1428 static void calc_column_indices(nbnxn_grid_t *grid,
1429                                 int a0, int a1,
1430                                 rvec *x,
1431                                 int dd_zone, const int *move,
1432                                 int thread, int nthread,
1433                                 int *cell,
1434                                 int *cxy_na)
1435 {
1436     int  n0, n1, i;
1437     int  cx, cy;
1438
1439     /* We add one extra cell for particles which moved during DD */
1440     for (i = 0; i < grid->ncx*grid->ncy+1; i++)
1441     {
1442         cxy_na[i] = 0;
1443     }
1444
1445     n0 = a0 + (int)((thread+0)*(a1 - a0))/nthread;
1446     n1 = a0 + (int)((thread+1)*(a1 - a0))/nthread;
1447     if (dd_zone == 0)
1448     {
1449         /* Home zone */
1450         for (i = n0; i < n1; i++)
1451         {
1452             if (move == NULL || move[i] >= 0)
1453             {
1454                 /* We need to be careful with rounding,
1455                  * particles might be a few bits outside the local zone.
1456                  * The int cast takes care of the lower bound,
1457                  * we will explicitly take care of the upper bound.
1458                  */
1459                 cx = (int)((x[i][XX] - grid->c0[XX])*grid->inv_sx);
1460                 cy = (int)((x[i][YY] - grid->c0[YY])*grid->inv_sy);
1461
1462 #ifndef NDEBUG
1463                 if (cx < 0 || cx > grid->ncx ||
1464                     cy < 0 || cy > grid->ncy)
1465                 {
1466                     gmx_fatal(FARGS,
1467                               "grid cell cx %d cy %d out of range (max %d %d)\n"
1468                               "atom %f %f %f, grid->c0 %f %f",
1469                               cx, cy, grid->ncx, grid->ncy,
1470                               x[i][XX], x[i][YY], x[i][ZZ], grid->c0[XX], grid->c0[YY]);
1471                 }
1472 #endif
1473                 /* Take care of potential rouding issues */
1474                 cx = min(cx, grid->ncx - 1);
1475                 cy = min(cy, grid->ncy - 1);
1476
1477                 /* For the moment cell will contain only the, grid local,
1478                  * x and y indices, not z.
1479                  */
1480                 cell[i] = cx*grid->ncy + cy;
1481             }
1482             else
1483             {
1484                 /* Put this moved particle after the end of the grid,
1485                  * so we can process it later without using conditionals.
1486                  */
1487                 cell[i] = grid->ncx*grid->ncy;
1488             }
1489
1490             cxy_na[cell[i]]++;
1491         }
1492     }
1493     else
1494     {
1495         /* Non-home zone */
1496         for (i = n0; i < n1; i++)
1497         {
1498             cx = (int)((x[i][XX] - grid->c0[XX])*grid->inv_sx);
1499             cy = (int)((x[i][YY] - grid->c0[YY])*grid->inv_sy);
1500
1501             /* For non-home zones there could be particles outside
1502              * the non-bonded cut-off range, which have been communicated
1503              * for bonded interactions only. For the result it doesn't
1504              * matter where these end up on the grid. For performance
1505              * we put them in an extra row at the border.
1506              */
1507             cx = max(cx, 0);
1508             cx = min(cx, grid->ncx - 1);
1509             cy = max(cy, 0);
1510             cy = min(cy, grid->ncy - 1);
1511
1512             /* For the moment cell will contain only the, grid local,
1513              * x and y indices, not z.
1514              */
1515             cell[i] = cx*grid->ncy + cy;
1516
1517             cxy_na[cell[i]]++;
1518         }
1519     }
1520 }
1521
1522 /* Determine in which grid cells the atoms should go */
1523 static void calc_cell_indices(const nbnxn_search_t nbs,
1524                               int dd_zone,
1525                               nbnxn_grid_t *grid,
1526                               int a0, int a1,
1527                               const int *atinfo,
1528                               rvec *x,
1529                               const int *move,
1530                               nbnxn_atomdata_t *nbat)
1531 {
1532     int   n0, n1, i;
1533     int   cx, cy, cxy, ncz_max, ncz;
1534     int   nthread, thread;
1535     int  *cxy_na, cxy_na_i;
1536
1537     nthread = gmx_omp_nthreads_get(emntPairsearch);
1538
1539 #pragma omp parallel for num_threads(nthread) schedule(static)
1540     for (thread = 0; thread < nthread; thread++)
1541     {
1542         calc_column_indices(grid, a0, a1, x, dd_zone, move, thread, nthread,
1543                             nbs->cell, nbs->work[thread].cxy_na);
1544     }
1545
1546     /* Make the cell index as a function of x and y */
1547     ncz_max          = 0;
1548     ncz              = 0;
1549     grid->cxy_ind[0] = 0;
1550     for (i = 0; i < grid->ncx*grid->ncy+1; i++)
1551     {
1552         /* We set ncz_max at the beginning of the loop iso at the end
1553          * to skip i=grid->ncx*grid->ncy which are moved particles
1554          * that do not need to be ordered on the grid.
1555          */
1556         if (ncz > ncz_max)
1557         {
1558             ncz_max = ncz;
1559         }
1560         cxy_na_i = nbs->work[0].cxy_na[i];
1561         for (thread = 1; thread < nthread; thread++)
1562         {
1563             cxy_na_i += nbs->work[thread].cxy_na[i];
1564         }
1565         ncz = (cxy_na_i + grid->na_sc - 1)/grid->na_sc;
1566         if (nbat->XFormat == nbatX8)
1567         {
1568             /* Make the number of cell a multiple of 2 */
1569             ncz = (ncz + 1) & ~1;
1570         }
1571         grid->cxy_ind[i+1] = grid->cxy_ind[i] + ncz;
1572         /* Clear cxy_na, so we can reuse the array below */
1573         grid->cxy_na[i] = 0;
1574     }
1575     grid->nc = grid->cxy_ind[grid->ncx*grid->ncy] - grid->cxy_ind[0];
1576
1577     nbat->natoms = (grid->cell0 + grid->nc)*grid->na_sc;
1578
1579     if (debug)
1580     {
1581         fprintf(debug, "ns na_sc %d na_c %d super-cells: %d x %d y %d z %.1f maxz %d\n",
1582                 grid->na_sc, grid->na_c, grid->nc,
1583                 grid->ncx, grid->ncy, grid->nc/((double)(grid->ncx*grid->ncy)),
1584                 ncz_max);
1585         if (gmx_debug_at)
1586         {
1587             i = 0;
1588             for (cy = 0; cy < grid->ncy; cy++)
1589             {
1590                 for (cx = 0; cx < grid->ncx; cx++)
1591                 {
1592                     fprintf(debug, " %2d", grid->cxy_ind[i+1]-grid->cxy_ind[i]);
1593                     i++;
1594                 }
1595                 fprintf(debug, "\n");
1596             }
1597         }
1598     }
1599
1600     /* Make sure the work array for sorting is large enough */
1601     if (ncz_max*grid->na_sc*SGSF > nbs->work[0].sort_work_nalloc)
1602     {
1603         for (thread = 0; thread < nbs->nthread_max; thread++)
1604         {
1605             nbs->work[thread].sort_work_nalloc =
1606                 over_alloc_large(ncz_max*grid->na_sc*SGSF);
1607             srenew(nbs->work[thread].sort_work,
1608                    nbs->work[thread].sort_work_nalloc);
1609             /* When not in use, all elements should be -1 */
1610             for (i = 0; i < nbs->work[thread].sort_work_nalloc; i++)
1611             {
1612                 nbs->work[thread].sort_work[i] = -1;
1613             }
1614         }
1615     }
1616
1617     /* Now we know the dimensions we can fill the grid.
1618      * This is the first, unsorted fill. We sort the columns after this.
1619      */
1620     for (i = a0; i < a1; i++)
1621     {
1622         /* At this point nbs->cell contains the local grid x,y indices */
1623         cxy = nbs->cell[i];
1624         nbs->a[(grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc + grid->cxy_na[cxy]++] = i;
1625     }
1626
1627     if (dd_zone == 0)
1628     {
1629         /* Set the cell indices for the moved particles */
1630         n0 = grid->nc*grid->na_sc;
1631         n1 = grid->nc*grid->na_sc+grid->cxy_na[grid->ncx*grid->ncy];
1632         if (dd_zone == 0)
1633         {
1634             for (i = n0; i < n1; i++)
1635             {
1636                 nbs->cell[nbs->a[i]] = i;
1637             }
1638         }
1639     }
1640
1641     /* Sort the super-cell columns along z into the sub-cells. */
1642 #pragma omp parallel for num_threads(nbs->nthread_max) schedule(static)
1643     for (thread = 0; thread < nbs->nthread_max; thread++)
1644     {
1645         if (grid->bSimple)
1646         {
1647             sort_columns_simple(nbs, dd_zone, grid, a0, a1, atinfo, x, nbat,
1648                                 ((thread+0)*grid->ncx*grid->ncy)/nthread,
1649                                 ((thread+1)*grid->ncx*grid->ncy)/nthread,
1650                                 nbs->work[thread].sort_work);
1651         }
1652         else
1653         {
1654             sort_columns_supersub(nbs, dd_zone, grid, a0, a1, atinfo, x, nbat,
1655                                   ((thread+0)*grid->ncx*grid->ncy)/nthread,
1656                                   ((thread+1)*grid->ncx*grid->ncy)/nthread,
1657                                   nbs->work[thread].sort_work);
1658         }
1659     }
1660
1661     if (grid->bSimple && nbat->XFormat == nbatX8)
1662     {
1663         combine_bounding_box_pairs(grid, grid->bb);
1664     }
1665
1666     if (!grid->bSimple)
1667     {
1668         grid->nsubc_tot = 0;
1669         for (i = 0; i < grid->nc; i++)
1670         {
1671             grid->nsubc_tot += grid->nsubc[i];
1672         }
1673     }
1674
1675     if (debug)
1676     {
1677         if (grid->bSimple)
1678         {
1679             print_bbsizes_simple(debug, nbs, grid);
1680         }
1681         else
1682         {
1683             fprintf(debug, "ns non-zero sub-cells: %d average atoms %.2f\n",
1684                     grid->nsubc_tot, (a1-a0)/(double)grid->nsubc_tot);
1685
1686             print_bbsizes_supersub(debug, nbs, grid);
1687         }
1688     }
1689 }
1690
1691 static void init_buffer_flags(nbnxn_buffer_flags_t *flags,
1692                               int                   natoms)
1693 {
1694     int b;
1695
1696     flags->nflag = (natoms + NBNXN_BUFFERFLAG_SIZE - 1)/NBNXN_BUFFERFLAG_SIZE;
1697     if (flags->nflag > flags->flag_nalloc)
1698     {
1699         flags->flag_nalloc = over_alloc_large(flags->nflag);
1700         srenew(flags->flag, flags->flag_nalloc);
1701     }
1702     for (b = 0; b < flags->nflag; b++)
1703     {
1704         flags->flag[b] = 0;
1705     }
1706 }
1707
1708 /* Sets up a grid and puts the atoms on the grid.
1709  * This function only operates on one domain of the domain decompostion.
1710  * Note that without domain decomposition there is only one domain.
1711  */
1712 void nbnxn_put_on_grid(nbnxn_search_t nbs,
1713                        int ePBC, matrix box,
1714                        int dd_zone,
1715                        rvec corner0, rvec corner1,
1716                        int a0, int a1,
1717                        real atom_density,
1718                        const int *atinfo,
1719                        rvec *x,
1720                        int nmoved, int *move,
1721                        int nb_kernel_type,
1722                        nbnxn_atomdata_t *nbat)
1723 {
1724     nbnxn_grid_t *grid;
1725     int           n;
1726     int           nc_max_grid, nc_max;
1727
1728     grid = &nbs->grid[dd_zone];
1729
1730     nbs_cycle_start(&nbs->cc[enbsCCgrid]);
1731
1732     grid->bSimple = nbnxn_kernel_pairlist_simple(nb_kernel_type);
1733
1734     grid->na_c      = nbnxn_kernel_to_ci_size(nb_kernel_type);
1735     grid->na_cj     = nbnxn_kernel_to_cj_size(nb_kernel_type);
1736     grid->na_sc     = (grid->bSimple ? 1 : GPU_NSUBCELL)*grid->na_c;
1737     grid->na_c_2log = get_2log(grid->na_c);
1738
1739     nbat->na_c = grid->na_c;
1740
1741     if (dd_zone == 0)
1742     {
1743         grid->cell0 = 0;
1744     }
1745     else
1746     {
1747         grid->cell0 =
1748             (nbs->grid[dd_zone-1].cell0 + nbs->grid[dd_zone-1].nc)*
1749             nbs->grid[dd_zone-1].na_sc/grid->na_sc;
1750     }
1751
1752     n = a1 - a0;
1753
1754     if (dd_zone == 0)
1755     {
1756         nbs->ePBC = ePBC;
1757         copy_mat(box, nbs->box);
1758
1759         if (atom_density >= 0)
1760         {
1761             grid->atom_density = atom_density;
1762         }
1763         else
1764         {
1765             grid->atom_density = grid_atom_density(n-nmoved, corner0, corner1);
1766         }
1767
1768         grid->cell0 = 0;
1769
1770         nbs->natoms_local    = a1 - nmoved;
1771         /* We assume that nbnxn_put_on_grid is called first
1772          * for the local atoms (dd_zone=0).
1773          */
1774         nbs->natoms_nonlocal = a1 - nmoved;
1775     }
1776     else
1777     {
1778         nbs->natoms_nonlocal = max(nbs->natoms_nonlocal, a1);
1779     }
1780
1781     nc_max_grid = set_grid_size_xy(nbs, grid,
1782                                    dd_zone, n-nmoved, corner0, corner1,
1783                                    nbs->grid[0].atom_density);
1784
1785     nc_max = grid->cell0 + nc_max_grid;
1786
1787     if (a1 > nbs->cell_nalloc)
1788     {
1789         nbs->cell_nalloc = over_alloc_large(a1);
1790         srenew(nbs->cell, nbs->cell_nalloc);
1791     }
1792
1793     /* To avoid conditionals we store the moved particles at the end of a,
1794      * make sure we have enough space.
1795      */
1796     if (nc_max*grid->na_sc + nmoved > nbs->a_nalloc)
1797     {
1798         nbs->a_nalloc = over_alloc_large(nc_max*grid->na_sc + nmoved);
1799         srenew(nbs->a, nbs->a_nalloc);
1800     }
1801
1802     /* We need padding up to a multiple of the buffer flag size: simply add */
1803     if (nc_max*grid->na_sc + NBNXN_BUFFERFLAG_SIZE > nbat->nalloc)
1804     {
1805         nbnxn_atomdata_realloc(nbat, nc_max*grid->na_sc+NBNXN_BUFFERFLAG_SIZE);
1806     }
1807
1808     calc_cell_indices(nbs, dd_zone, grid, a0, a1, atinfo, x, move, nbat);
1809
1810     if (dd_zone == 0)
1811     {
1812         nbat->natoms_local = nbat->natoms;
1813     }
1814
1815     nbs_cycle_stop(&nbs->cc[enbsCCgrid]);
1816 }
1817
1818 /* Calls nbnxn_put_on_grid for all non-local domains */
1819 void nbnxn_put_on_grid_nonlocal(nbnxn_search_t            nbs,
1820                                 const gmx_domdec_zones_t *zones,
1821                                 const int                *atinfo,
1822                                 rvec                     *x,
1823                                 int                       nb_kernel_type,
1824                                 nbnxn_atomdata_t         *nbat)
1825 {
1826     int  zone, d;
1827     rvec c0, c1;
1828
1829     for (zone = 1; zone < zones->n; zone++)
1830     {
1831         for (d = 0; d < DIM; d++)
1832         {
1833             c0[d] = zones->size[zone].bb_x0[d];
1834             c1[d] = zones->size[zone].bb_x1[d];
1835         }
1836
1837         nbnxn_put_on_grid(nbs, nbs->ePBC, NULL,
1838                           zone, c0, c1,
1839                           zones->cg_range[zone],
1840                           zones->cg_range[zone+1],
1841                           -1,
1842                           atinfo,
1843                           x,
1844                           0, NULL,
1845                           nb_kernel_type,
1846                           nbat);
1847     }
1848 }
1849
1850 /* Add simple grid type information to the local super/sub grid */
1851 void nbnxn_grid_add_simple(nbnxn_search_t    nbs,
1852                            nbnxn_atomdata_t *nbat)
1853 {
1854     nbnxn_grid_t *grid;
1855     float        *bbcz;
1856     nbnxn_bb_t   *bb;
1857     int           ncd, sc;
1858
1859     grid = &nbs->grid[0];
1860
1861     if (grid->bSimple)
1862     {
1863         gmx_incons("nbnxn_grid_simple called with a simple grid");
1864     }
1865
1866     ncd = grid->na_sc/NBNXN_CPU_CLUSTER_I_SIZE;
1867
1868     if (grid->nc*ncd > grid->nc_nalloc_simple)
1869     {
1870         grid->nc_nalloc_simple = over_alloc_large(grid->nc*ncd);
1871         srenew(grid->bbcz_simple, grid->nc_nalloc_simple*NNBSBB_D);
1872         srenew(grid->bb_simple, grid->nc_nalloc_simple);
1873         srenew(grid->flags_simple, grid->nc_nalloc_simple);
1874         if (nbat->XFormat)
1875         {
1876             sfree_aligned(grid->bbj);
1877             snew_aligned(grid->bbj, grid->nc_nalloc_simple/2, 16);
1878         }
1879     }
1880
1881     bbcz = grid->bbcz_simple;
1882     bb   = grid->bb_simple;
1883
1884 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntPairsearch)) schedule(static)
1885     for (sc = 0; sc < grid->nc; sc++)
1886     {
1887         int c, tx, na;
1888
1889         for (c = 0; c < ncd; c++)
1890         {
1891             tx = sc*ncd + c;
1892
1893             na = NBNXN_CPU_CLUSTER_I_SIZE;
1894             while (na > 0 &&
1895                    nbat->type[tx*NBNXN_CPU_CLUSTER_I_SIZE+na-1] == nbat->ntype-1)
1896             {
1897                 na--;
1898             }
1899
1900             if (na > 0)
1901             {
1902                 switch (nbat->XFormat)
1903                 {
1904                     case nbatX4:
1905                         /* PACK_X4==NBNXN_CPU_CLUSTER_I_SIZE, so this is simple */
1906                         calc_bounding_box_x_x4(na, nbat->x+tx*STRIDE_P4,
1907                                                bb+tx);
1908                         break;
1909                     case nbatX8:
1910                         /* PACK_X8>NBNXN_CPU_CLUSTER_I_SIZE, more complicated */
1911                         calc_bounding_box_x_x8(na, nbat->x+X8_IND_A(tx*NBNXN_CPU_CLUSTER_I_SIZE),
1912                                                bb+tx);
1913                         break;
1914                     default:
1915                         calc_bounding_box(na, nbat->xstride,
1916                                           nbat->x+tx*NBNXN_CPU_CLUSTER_I_SIZE*nbat->xstride,
1917                                           bb+tx);
1918                         break;
1919                 }
1920                 bbcz[tx*NNBSBB_D+0] = bb[tx].lower[BB_Z];
1921                 bbcz[tx*NNBSBB_D+1] = bb[tx].upper[BB_Z];
1922
1923                 /* No interaction optimization yet here */
1924                 grid->flags_simple[tx] = NBNXN_CI_DO_LJ(0) | NBNXN_CI_DO_COUL(0);
1925             }
1926             else
1927             {
1928                 grid->flags_simple[tx] = 0;
1929             }
1930         }
1931     }
1932
1933     if (grid->bSimple && nbat->XFormat == nbatX8)
1934     {
1935         combine_bounding_box_pairs(grid, grid->bb_simple);
1936     }
1937 }
1938
1939 void nbnxn_get_ncells(nbnxn_search_t nbs, int *ncx, int *ncy)
1940 {
1941     *ncx = nbs->grid[0].ncx;
1942     *ncy = nbs->grid[0].ncy;
1943 }
1944
1945 void nbnxn_get_atomorder(nbnxn_search_t nbs, int **a, int *n)
1946 {
1947     const nbnxn_grid_t *grid;
1948
1949     grid = &nbs->grid[0];
1950
1951     /* Return the atom order for the home cell (index 0) */
1952     *a  = nbs->a;
1953
1954     *n = grid->cxy_ind[grid->ncx*grid->ncy]*grid->na_sc;
1955 }
1956
1957 void nbnxn_set_atomorder(nbnxn_search_t nbs)
1958 {
1959     nbnxn_grid_t *grid;
1960     int           ao, cx, cy, cxy, cz, j;
1961
1962     /* Set the atom order for the home cell (index 0) */
1963     grid = &nbs->grid[0];
1964
1965     ao = 0;
1966     for (cx = 0; cx < grid->ncx; cx++)
1967     {
1968         for (cy = 0; cy < grid->ncy; cy++)
1969         {
1970             cxy = cx*grid->ncy + cy;
1971             j   = grid->cxy_ind[cxy]*grid->na_sc;
1972             for (cz = 0; cz < grid->cxy_na[cxy]; cz++)
1973             {
1974                 nbs->a[j]     = ao;
1975                 nbs->cell[ao] = j;
1976                 ao++;
1977                 j++;
1978             }
1979         }
1980     }
1981 }
1982
1983 /* Determines the cell range along one dimension that
1984  * the bounding box b0 - b1 sees.
1985  */
1986 static void get_cell_range(real b0, real b1,
1987                            int nc, real c0, real s, real invs,
1988                            real d2, real r2, int *cf, int *cl)
1989 {
1990     *cf = max((int)((b0 - c0)*invs), 0);
1991
1992     while (*cf > 0 && d2 + sqr((b0 - c0) - (*cf-1+1)*s) < r2)
1993     {
1994         (*cf)--;
1995     }
1996
1997     *cl = min((int)((b1 - c0)*invs), nc-1);
1998     while (*cl < nc-1 && d2 + sqr((*cl+1)*s - (b1 - c0)) < r2)
1999     {
2000         (*cl)++;
2001     }
2002 }
2003
2004 /* Reference code calculating the distance^2 between two bounding boxes */
2005 static float box_dist2(float bx0, float bx1, float by0,
2006                        float by1, float bz0, float bz1,
2007                        const nbnxn_bb_t *bb)
2008 {
2009     float d2;
2010     float dl, dh, dm, dm0;
2011
2012     d2 = 0;
2013
2014     dl  = bx0 - bb->upper[BB_X];
2015     dh  = bb->lower[BB_X] - bx1;
2016     dm  = max(dl, dh);
2017     dm0 = max(dm, 0);
2018     d2 += dm0*dm0;
2019
2020     dl  = by0 - bb->upper[BB_Y];
2021     dh  = bb->lower[BB_Y] - by1;
2022     dm  = max(dl, dh);
2023     dm0 = max(dm, 0);
2024     d2 += dm0*dm0;
2025
2026     dl  = bz0 - bb->upper[BB_Z];
2027     dh  = bb->lower[BB_Z] - bz1;
2028     dm  = max(dl, dh);
2029     dm0 = max(dm, 0);
2030     d2 += dm0*dm0;
2031
2032     return d2;
2033 }
2034
2035 /* Plain C code calculating the distance^2 between two bounding boxes */
2036 static float subc_bb_dist2(int si, const nbnxn_bb_t *bb_i_ci,
2037                            int csj, const nbnxn_bb_t *bb_j_all)
2038 {
2039     const nbnxn_bb_t *bb_i, *bb_j;
2040     float             d2;
2041     float             dl, dh, dm, dm0;
2042
2043     bb_i = bb_i_ci  +  si;
2044     bb_j = bb_j_all + csj;
2045
2046     d2 = 0;
2047
2048     dl  = bb_i->lower[BB_X] - bb_j->upper[BB_X];
2049     dh  = bb_j->lower[BB_X] - bb_i->upper[BB_X];
2050     dm  = max(dl, dh);
2051     dm0 = max(dm, 0);
2052     d2 += dm0*dm0;
2053
2054     dl  = bb_i->lower[BB_Y] - bb_j->upper[BB_Y];
2055     dh  = bb_j->lower[BB_Y] - bb_i->upper[BB_Y];
2056     dm  = max(dl, dh);
2057     dm0 = max(dm, 0);
2058     d2 += dm0*dm0;
2059
2060     dl  = bb_i->lower[BB_Z] - bb_j->upper[BB_Z];
2061     dh  = bb_j->lower[BB_Z] - bb_i->upper[BB_Z];
2062     dm  = max(dl, dh);
2063     dm0 = max(dm, 0);
2064     d2 += dm0*dm0;
2065
2066     return d2;
2067 }
2068
2069 #ifdef NBNXN_SEARCH_BB_SIMD4
2070
2071 /* 4-wide SIMD code for bb distance for bb format xyz0 */
2072 static float subc_bb_dist2_simd4(int si, const nbnxn_bb_t *bb_i_ci,
2073                                  int csj, const nbnxn_bb_t *bb_j_all)
2074 {
2075     gmx_simd4_pr bb_i_S0, bb_i_S1;
2076     gmx_simd4_pr bb_j_S0, bb_j_S1;
2077     gmx_simd4_pr dl_S;
2078     gmx_simd4_pr dh_S;
2079     gmx_simd4_pr dm_S;
2080     gmx_simd4_pr dm0_S;
2081
2082     bb_i_S0 = gmx_simd4_load_bb_pr(&bb_i_ci[si].lower[0]);
2083     bb_i_S1 = gmx_simd4_load_bb_pr(&bb_i_ci[si].upper[0]);
2084     bb_j_S0 = gmx_simd4_load_bb_pr(&bb_j_all[csj].lower[0]);
2085     bb_j_S1 = gmx_simd4_load_bb_pr(&bb_j_all[csj].upper[0]);
2086
2087     dl_S    = gmx_simd4_sub_pr(bb_i_S0, bb_j_S1);
2088     dh_S    = gmx_simd4_sub_pr(bb_j_S0, bb_i_S1);
2089
2090     dm_S    = gmx_simd4_max_pr(dl_S, dh_S);
2091     dm0_S   = gmx_simd4_max_pr(dm_S, gmx_simd4_setzero_pr());
2092
2093     return gmx_simd4_dotproduct3(dm0_S, dm0_S);
2094 }
2095
2096 /* Calculate bb bounding distances of bb_i[si,...,si+3] and store them in d2 */
2097 #define SUBC_BB_DIST2_SIMD4_XXXX_INNER(si, bb_i, d2) \
2098     {                                                \
2099         int          shi;                                  \
2100                                                  \
2101         gmx_simd4_pr dx_0, dy_0, dz_0;                       \
2102         gmx_simd4_pr dx_1, dy_1, dz_1;                       \
2103                                                  \
2104         gmx_simd4_pr mx, my, mz;                             \
2105         gmx_simd4_pr m0x, m0y, m0z;                          \
2106                                                  \
2107         gmx_simd4_pr d2x, d2y, d2z;                          \
2108         gmx_simd4_pr d2s, d2t;                              \
2109                                                  \
2110         shi = si*NNBSBB_D*DIM;                       \
2111                                                  \
2112         xi_l = gmx_simd4_load_bb_pr(bb_i+shi+0*STRIDE_PBB);   \
2113         yi_l = gmx_simd4_load_bb_pr(bb_i+shi+1*STRIDE_PBB);   \
2114         zi_l = gmx_simd4_load_bb_pr(bb_i+shi+2*STRIDE_PBB);   \
2115         xi_h = gmx_simd4_load_bb_pr(bb_i+shi+3*STRIDE_PBB);   \
2116         yi_h = gmx_simd4_load_bb_pr(bb_i+shi+4*STRIDE_PBB);   \
2117         zi_h = gmx_simd4_load_bb_pr(bb_i+shi+5*STRIDE_PBB);   \
2118                                                  \
2119         dx_0 = gmx_simd4_sub_pr(xi_l, xj_h);                \
2120         dy_0 = gmx_simd4_sub_pr(yi_l, yj_h);                \
2121         dz_0 = gmx_simd4_sub_pr(zi_l, zj_h);                \
2122                                                  \
2123         dx_1 = gmx_simd4_sub_pr(xj_l, xi_h);                \
2124         dy_1 = gmx_simd4_sub_pr(yj_l, yi_h);                \
2125         dz_1 = gmx_simd4_sub_pr(zj_l, zi_h);                \
2126                                                  \
2127         mx   = gmx_simd4_max_pr(dx_0, dx_1);                \
2128         my   = gmx_simd4_max_pr(dy_0, dy_1);                \
2129         mz   = gmx_simd4_max_pr(dz_0, dz_1);                \
2130                                                  \
2131         m0x  = gmx_simd4_max_pr(mx, zero);                  \
2132         m0y  = gmx_simd4_max_pr(my, zero);                  \
2133         m0z  = gmx_simd4_max_pr(mz, zero);                  \
2134                                                  \
2135         d2x  = gmx_simd4_mul_pr(m0x, m0x);                  \
2136         d2y  = gmx_simd4_mul_pr(m0y, m0y);                  \
2137         d2z  = gmx_simd4_mul_pr(m0z, m0z);                  \
2138                                                  \
2139         d2s  = gmx_simd4_add_pr(d2x, d2y);                  \
2140         d2t  = gmx_simd4_add_pr(d2s, d2z);                  \
2141                                                  \
2142         gmx_simd4_store_pr(d2+si, d2t);                     \
2143     }
2144
2145 /* 4-wide SIMD code for nsi bb distances for bb format xxxxyyyyzzzz */
2146 static void subc_bb_dist2_simd4_xxxx(const float *bb_j,
2147                                      int nsi, const float *bb_i,
2148                                      float *d2)
2149 {
2150     gmx_simd4_pr xj_l, yj_l, zj_l;
2151     gmx_simd4_pr xj_h, yj_h, zj_h;
2152     gmx_simd4_pr xi_l, yi_l, zi_l;
2153     gmx_simd4_pr xi_h, yi_h, zi_h;
2154
2155     gmx_simd4_pr zero;
2156
2157     zero = gmx_simd4_setzero_pr();
2158
2159     xj_l = gmx_simd4_set1_pr(bb_j[0*STRIDE_PBB]);
2160     yj_l = gmx_simd4_set1_pr(bb_j[1*STRIDE_PBB]);
2161     zj_l = gmx_simd4_set1_pr(bb_j[2*STRIDE_PBB]);
2162     xj_h = gmx_simd4_set1_pr(bb_j[3*STRIDE_PBB]);
2163     yj_h = gmx_simd4_set1_pr(bb_j[4*STRIDE_PBB]);
2164     zj_h = gmx_simd4_set1_pr(bb_j[5*STRIDE_PBB]);
2165
2166     /* Here we "loop" over si (0,STRIDE_PBB) from 0 to nsi with step STRIDE_PBB.
2167      * But as we know the number of iterations is 1 or 2, we unroll manually.
2168      */
2169     SUBC_BB_DIST2_SIMD4_XXXX_INNER(0, bb_i, d2);
2170     if (STRIDE_PBB < nsi)
2171     {
2172         SUBC_BB_DIST2_SIMD4_XXXX_INNER(STRIDE_PBB, bb_i, d2);
2173     }
2174 }
2175
2176 #endif /* NBNXN_SEARCH_BB_SIMD4 */
2177
2178 /* Plain C function which determines if any atom pair between two cells
2179  * is within distance sqrt(rl2).
2180  */
2181 static gmx_bool subc_in_range_x(int na_c,
2182                                 int si, const real *x_i,
2183                                 int csj, int stride, const real *x_j,
2184                                 real rl2)
2185 {
2186     int  i, j, i0, j0;
2187     real d2;
2188
2189     for (i = 0; i < na_c; i++)
2190     {
2191         i0 = (si*na_c + i)*DIM;
2192         for (j = 0; j < na_c; j++)
2193         {
2194             j0 = (csj*na_c + j)*stride;
2195
2196             d2 = sqr(x_i[i0  ] - x_j[j0  ]) +
2197                 sqr(x_i[i0+1] - x_j[j0+1]) +
2198                 sqr(x_i[i0+2] - x_j[j0+2]);
2199
2200             if (d2 < rl2)
2201             {
2202                 return TRUE;
2203             }
2204         }
2205     }
2206
2207     return FALSE;
2208 }
2209
2210 #ifdef NBNXN_SEARCH_SIMD4_FLOAT_X_BB
2211 /* When we make seperate single/double precision SIMD vector operation
2212  * include files, this function should be moved there (also using FMA).
2213  */
2214 static inline gmx_simd4_pr
2215 gmx_simd4_calc_rsq_pr(gmx_simd4_pr x, gmx_simd4_pr y, gmx_simd4_pr z)
2216 {
2217     return gmx_simd4_add_pr( gmx_simd4_add_pr( gmx_simd4_mul_pr(x, x), gmx_simd4_mul_pr(y, y) ), gmx_simd4_mul_pr(z, z) );
2218 }
2219
2220 /* 4-wide SIMD function which determines if any atom pair between two cells,
2221  * both with 8 atoms, is within distance sqrt(rl2).
2222  * Using 8-wide AVX is not faster on Intel Sandy Bridge.
2223  */
2224 static gmx_bool subc_in_range_simd4(int na_c,
2225                                     int si, const real *x_i,
2226                                     int csj, int stride, const real *x_j,
2227                                     real rl2)
2228 {
2229     gmx_simd4_pr ix_S0, iy_S0, iz_S0;
2230     gmx_simd4_pr ix_S1, iy_S1, iz_S1;
2231
2232     gmx_simd4_pr rc2_S;
2233
2234     int          dim_stride;
2235     int          j0, j1;
2236
2237     rc2_S   = gmx_simd4_set1_pr(rl2);
2238
2239     dim_stride = NBNXN_GPU_CLUSTER_SIZE/STRIDE_PBB*DIM;
2240     ix_S0      = gmx_simd4_load_bb_pr(x_i+(si*dim_stride+0)*STRIDE_PBB);
2241     iy_S0      = gmx_simd4_load_bb_pr(x_i+(si*dim_stride+1)*STRIDE_PBB);
2242     iz_S0      = gmx_simd4_load_bb_pr(x_i+(si*dim_stride+2)*STRIDE_PBB);
2243     ix_S1      = gmx_simd4_load_bb_pr(x_i+(si*dim_stride+3)*STRIDE_PBB);
2244     iy_S1      = gmx_simd4_load_bb_pr(x_i+(si*dim_stride+4)*STRIDE_PBB);
2245     iz_S1      = gmx_simd4_load_bb_pr(x_i+(si*dim_stride+5)*STRIDE_PBB);
2246
2247     /* We loop from the outer to the inner particles to maximize
2248      * the chance that we find a pair in range quickly and return.
2249      */
2250     j0 = csj*na_c;
2251     j1 = j0 + na_c - 1;
2252     while (j0 < j1)
2253     {
2254         gmx_simd4_pr jx0_S, jy0_S, jz0_S;
2255         gmx_simd4_pr jx1_S, jy1_S, jz1_S;
2256
2257         gmx_simd4_pr dx_S0, dy_S0, dz_S0;
2258         gmx_simd4_pr dx_S1, dy_S1, dz_S1;
2259         gmx_simd4_pr dx_S2, dy_S2, dz_S2;
2260         gmx_simd4_pr dx_S3, dy_S3, dz_S3;
2261
2262         gmx_simd4_pr rsq_S0;
2263         gmx_simd4_pr rsq_S1;
2264         gmx_simd4_pr rsq_S2;
2265         gmx_simd4_pr rsq_S3;
2266
2267         gmx_simd4_pb wco_S0;
2268         gmx_simd4_pb wco_S1;
2269         gmx_simd4_pb wco_S2;
2270         gmx_simd4_pb wco_S3;
2271         gmx_simd4_pb wco_any_S01, wco_any_S23, wco_any_S;
2272
2273         jx0_S = gmx_simd4_set1_pr(x_j[j0*stride+0]);
2274         jy0_S = gmx_simd4_set1_pr(x_j[j0*stride+1]);
2275         jz0_S = gmx_simd4_set1_pr(x_j[j0*stride+2]);
2276
2277         jx1_S = gmx_simd4_set1_pr(x_j[j1*stride+0]);
2278         jy1_S = gmx_simd4_set1_pr(x_j[j1*stride+1]);
2279         jz1_S = gmx_simd4_set1_pr(x_j[j1*stride+2]);
2280
2281         /* Calculate distance */
2282         dx_S0            = gmx_simd4_sub_pr(ix_S0, jx0_S);
2283         dy_S0            = gmx_simd4_sub_pr(iy_S0, jy0_S);
2284         dz_S0            = gmx_simd4_sub_pr(iz_S0, jz0_S);
2285         dx_S1            = gmx_simd4_sub_pr(ix_S1, jx0_S);
2286         dy_S1            = gmx_simd4_sub_pr(iy_S1, jy0_S);
2287         dz_S1            = gmx_simd4_sub_pr(iz_S1, jz0_S);
2288         dx_S2            = gmx_simd4_sub_pr(ix_S0, jx1_S);
2289         dy_S2            = gmx_simd4_sub_pr(iy_S0, jy1_S);
2290         dz_S2            = gmx_simd4_sub_pr(iz_S0, jz1_S);
2291         dx_S3            = gmx_simd4_sub_pr(ix_S1, jx1_S);
2292         dy_S3            = gmx_simd4_sub_pr(iy_S1, jy1_S);
2293         dz_S3            = gmx_simd4_sub_pr(iz_S1, jz1_S);
2294
2295         /* rsq = dx*dx+dy*dy+dz*dz */
2296         rsq_S0           = gmx_simd4_calc_rsq_pr(dx_S0, dy_S0, dz_S0);
2297         rsq_S1           = gmx_simd4_calc_rsq_pr(dx_S1, dy_S1, dz_S1);
2298         rsq_S2           = gmx_simd4_calc_rsq_pr(dx_S2, dy_S2, dz_S2);
2299         rsq_S3           = gmx_simd4_calc_rsq_pr(dx_S3, dy_S3, dz_S3);
2300
2301         wco_S0           = gmx_simd4_cmplt_pr(rsq_S0, rc2_S);
2302         wco_S1           = gmx_simd4_cmplt_pr(rsq_S1, rc2_S);
2303         wco_S2           = gmx_simd4_cmplt_pr(rsq_S2, rc2_S);
2304         wco_S3           = gmx_simd4_cmplt_pr(rsq_S3, rc2_S);
2305
2306         wco_any_S01      = gmx_simd4_or_pb(wco_S0, wco_S1);
2307         wco_any_S23      = gmx_simd4_or_pb(wco_S2, wco_S3);
2308         wco_any_S        = gmx_simd4_or_pb(wco_any_S01, wco_any_S23);
2309
2310         if (gmx_simd4_anytrue_pb(wco_any_S))
2311         {
2312             return TRUE;
2313         }
2314
2315         j0++;
2316         j1--;
2317     }
2318     return FALSE;
2319
2320 }
2321 #endif
2322
2323
2324 /* Returns the j sub-cell for index cj_ind */
2325 static int nbl_cj(const nbnxn_pairlist_t *nbl, int cj_ind)
2326 {
2327     return nbl->cj4[cj_ind >> NBNXN_GPU_JGROUP_SIZE_2LOG].cj[cj_ind & (NBNXN_GPU_JGROUP_SIZE - 1)];
2328 }
2329
2330 /* Returns the i-interaction mask of the j sub-cell for index cj_ind */
2331 static unsigned nbl_imask0(const nbnxn_pairlist_t *nbl, int cj_ind)
2332 {
2333     return nbl->cj4[cj_ind >> NBNXN_GPU_JGROUP_SIZE_2LOG].imei[0].imask;
2334 }
2335
2336 /* Ensures there is enough space for extra extra exclusion masks */
2337 static void check_excl_space(nbnxn_pairlist_t *nbl, int extra)
2338 {
2339     if (nbl->nexcl+extra > nbl->excl_nalloc)
2340     {
2341         nbl->excl_nalloc = over_alloc_small(nbl->nexcl+extra);
2342         nbnxn_realloc_void((void **)&nbl->excl,
2343                            nbl->nexcl*sizeof(*nbl->excl),
2344                            nbl->excl_nalloc*sizeof(*nbl->excl),
2345                            nbl->alloc, nbl->free);
2346     }
2347 }
2348
2349 /* Ensures there is enough space for ncell extra j-cells in the list */
2350 static void check_subcell_list_space_simple(nbnxn_pairlist_t *nbl,
2351                                             int               ncell)
2352 {
2353     int cj_max;
2354
2355     cj_max = nbl->ncj + ncell;
2356
2357     if (cj_max > nbl->cj_nalloc)
2358     {
2359         nbl->cj_nalloc = over_alloc_small(cj_max);
2360         nbnxn_realloc_void((void **)&nbl->cj,
2361                            nbl->ncj*sizeof(*nbl->cj),
2362                            nbl->cj_nalloc*sizeof(*nbl->cj),
2363                            nbl->alloc, nbl->free);
2364     }
2365 }
2366
2367 /* Ensures there is enough space for ncell extra j-subcells in the list */
2368 static void check_subcell_list_space_supersub(nbnxn_pairlist_t *nbl,
2369                                               int               nsupercell)
2370 {
2371     int ncj4_max, j4, j, w, t;
2372
2373 #define NWARP       2
2374 #define WARP_SIZE  32
2375
2376     /* We can have maximally nsupercell*GPU_NSUBCELL sj lists */
2377     /* We can store 4 j-subcell - i-supercell pairs in one struct.
2378      * since we round down, we need one extra entry.
2379      */
2380     ncj4_max = ((nbl->work->cj_ind + nsupercell*GPU_NSUBCELL + NBNXN_GPU_JGROUP_SIZE - 1) >> NBNXN_GPU_JGROUP_SIZE_2LOG);
2381
2382     if (ncj4_max > nbl->cj4_nalloc)
2383     {
2384         nbl->cj4_nalloc = over_alloc_small(ncj4_max);
2385         nbnxn_realloc_void((void **)&nbl->cj4,
2386                            nbl->work->cj4_init*sizeof(*nbl->cj4),
2387                            nbl->cj4_nalloc*sizeof(*nbl->cj4),
2388                            nbl->alloc, nbl->free);
2389     }
2390
2391     if (ncj4_max > nbl->work->cj4_init)
2392     {
2393         for (j4 = nbl->work->cj4_init; j4 < ncj4_max; j4++)
2394         {
2395             /* No i-subcells and no excl's in the list initially */
2396             for (w = 0; w < NWARP; w++)
2397             {
2398                 nbl->cj4[j4].imei[w].imask    = 0U;
2399                 nbl->cj4[j4].imei[w].excl_ind = 0;
2400
2401             }
2402         }
2403         nbl->work->cj4_init = ncj4_max;
2404     }
2405 }
2406
2407 /* Set all excl masks for one GPU warp no exclusions */
2408 static void set_no_excls(nbnxn_excl_t *excl)
2409 {
2410     int t;
2411
2412     for (t = 0; t < WARP_SIZE; t++)
2413     {
2414         /* Turn all interaction bits on */
2415         excl->pair[t] = NBNXN_INTERACTION_MASK_ALL;
2416     }
2417 }
2418
2419 /* Initializes a single nbnxn_pairlist_t data structure */
2420 static void nbnxn_init_pairlist(nbnxn_pairlist_t *nbl,
2421                                 gmx_bool          bSimple,
2422                                 nbnxn_alloc_t    *alloc,
2423                                 nbnxn_free_t     *free)
2424 {
2425     if (alloc == NULL)
2426     {
2427         nbl->alloc = nbnxn_alloc_aligned;
2428     }
2429     else
2430     {
2431         nbl->alloc = alloc;
2432     }
2433     if (free == NULL)
2434     {
2435         nbl->free = nbnxn_free_aligned;
2436     }
2437     else
2438     {
2439         nbl->free = free;
2440     }
2441
2442     nbl->bSimple     = bSimple;
2443     nbl->na_sc       = 0;
2444     nbl->na_ci       = 0;
2445     nbl->na_cj       = 0;
2446     nbl->nci         = 0;
2447     nbl->ci          = NULL;
2448     nbl->ci_nalloc   = 0;
2449     nbl->ncj         = 0;
2450     nbl->cj          = NULL;
2451     nbl->cj_nalloc   = 0;
2452     nbl->ncj4        = 0;
2453     /* We need one element extra in sj, so alloc initially with 1 */
2454     nbl->cj4_nalloc  = 0;
2455     nbl->cj4         = NULL;
2456     nbl->nci_tot     = 0;
2457
2458     if (!nbl->bSimple)
2459     {
2460         nbl->excl        = NULL;
2461         nbl->excl_nalloc = 0;
2462         nbl->nexcl       = 0;
2463         check_excl_space(nbl, 1);
2464         nbl->nexcl       = 1;
2465         set_no_excls(&nbl->excl[0]);
2466     }
2467
2468     snew(nbl->work, 1);
2469     if (nbl->bSimple)
2470     {
2471         snew_aligned(nbl->work->bb_ci, 1, NBNXN_SEARCH_BB_MEM_ALIGN);
2472     }
2473     else
2474     {
2475 #ifdef NBNXN_BBXXXX
2476         snew_aligned(nbl->work->pbb_ci, GPU_NSUBCELL/STRIDE_PBB*NNBSBB_XXXX, NBNXN_SEARCH_BB_MEM_ALIGN);
2477 #else
2478         snew_aligned(nbl->work->bb_ci, GPU_NSUBCELL, NBNXN_SEARCH_BB_MEM_ALIGN);
2479 #endif
2480     }
2481     snew_aligned(nbl->work->x_ci, NBNXN_NA_SC_MAX*DIM, NBNXN_SEARCH_BB_MEM_ALIGN);
2482 #ifdef GMX_NBNXN_SIMD
2483     snew_aligned(nbl->work->x_ci_simd_4xn, 1, NBNXN_MEM_ALIGN);
2484     snew_aligned(nbl->work->x_ci_simd_2xnn, 1, NBNXN_MEM_ALIGN);
2485 #endif
2486     snew_aligned(nbl->work->d2, GPU_NSUBCELL, NBNXN_SEARCH_BB_MEM_ALIGN);
2487
2488     nbl->work->sort            = NULL;
2489     nbl->work->sort_nalloc     = 0;
2490     nbl->work->sci_sort        = NULL;
2491     nbl->work->sci_sort_nalloc = 0;
2492 }
2493
2494 void nbnxn_init_pairlist_set(nbnxn_pairlist_set_t *nbl_list,
2495                              gmx_bool bSimple, gmx_bool bCombined,
2496                              nbnxn_alloc_t *alloc,
2497                              nbnxn_free_t  *free)
2498 {
2499     int i;
2500
2501     nbl_list->bSimple   = bSimple;
2502     nbl_list->bCombined = bCombined;
2503
2504     nbl_list->nnbl = gmx_omp_nthreads_get(emntNonbonded);
2505
2506     if (!nbl_list->bCombined &&
2507         nbl_list->nnbl > NBNXN_BUFFERFLAG_MAX_THREADS)
2508     {
2509         gmx_fatal(FARGS, "%d OpenMP threads were requested. Since the non-bonded force buffer reduction is prohibitively slow with more than %d threads, we do not allow this. Use %d or less OpenMP threads.",
2510                   nbl_list->nnbl, NBNXN_BUFFERFLAG_MAX_THREADS, NBNXN_BUFFERFLAG_MAX_THREADS);
2511     }
2512
2513     snew(nbl_list->nbl, nbl_list->nnbl);
2514     /* Execute in order to avoid memory interleaving between threads */
2515 #pragma omp parallel for num_threads(nbl_list->nnbl) schedule(static)
2516     for (i = 0; i < nbl_list->nnbl; i++)
2517     {
2518         /* Allocate the nblist data structure locally on each thread
2519          * to optimize memory access for NUMA architectures.
2520          */
2521         snew(nbl_list->nbl[i], 1);
2522
2523         /* Only list 0 is used on the GPU, use normal allocation for i>0 */
2524         if (i == 0)
2525         {
2526             nbnxn_init_pairlist(nbl_list->nbl[i], nbl_list->bSimple, alloc, free);
2527         }
2528         else
2529         {
2530             nbnxn_init_pairlist(nbl_list->nbl[i], nbl_list->bSimple, NULL, NULL);
2531         }
2532     }
2533 }
2534
2535 /* Print statistics of a pair list, used for debug output */
2536 static void print_nblist_statistics_simple(FILE *fp, const nbnxn_pairlist_t *nbl,
2537                                            const nbnxn_search_t nbs, real rl)
2538 {
2539     const nbnxn_grid_t *grid;
2540     int                 cs[SHIFTS];
2541     int                 s, i, j;
2542     int                 npexcl;
2543
2544     /* This code only produces correct statistics with domain decomposition */
2545     grid = &nbs->grid[0];
2546
2547     fprintf(fp, "nbl nci %d ncj %d\n",
2548             nbl->nci, nbl->ncj);
2549     fprintf(fp, "nbl na_sc %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
2550             nbl->na_sc, rl, nbl->ncj, nbl->ncj/(double)grid->nc,
2551             nbl->ncj/(double)grid->nc*grid->na_sc,
2552             nbl->ncj/(double)grid->nc*grid->na_sc/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nc*grid->na_sc/det(nbs->box)));
2553
2554     fprintf(fp, "nbl average j cell list length %.1f\n",
2555             0.25*nbl->ncj/(double)nbl->nci);
2556
2557     for (s = 0; s < SHIFTS; s++)
2558     {
2559         cs[s] = 0;
2560     }
2561     npexcl = 0;
2562     for (i = 0; i < nbl->nci; i++)
2563     {
2564         cs[nbl->ci[i].shift & NBNXN_CI_SHIFT] +=
2565             nbl->ci[i].cj_ind_end - nbl->ci[i].cj_ind_start;
2566
2567         j = nbl->ci[i].cj_ind_start;
2568         while (j < nbl->ci[i].cj_ind_end &&
2569                nbl->cj[j].excl != NBNXN_INTERACTION_MASK_ALL)
2570         {
2571             npexcl++;
2572             j++;
2573         }
2574     }
2575     fprintf(fp, "nbl cell pairs, total: %d excl: %d %.1f%%\n",
2576             nbl->ncj, npexcl, 100*npexcl/(double)nbl->ncj);
2577     for (s = 0; s < SHIFTS; s++)
2578     {
2579         if (cs[s] > 0)
2580         {
2581             fprintf(fp, "nbl shift %2d ncj %3d\n", s, cs[s]);
2582         }
2583     }
2584 }
2585
2586 /* Print statistics of a pair lists, used for debug output */
2587 static void print_nblist_statistics_supersub(FILE *fp, const nbnxn_pairlist_t *nbl,
2588                                              const nbnxn_search_t nbs, real rl)
2589 {
2590     const nbnxn_grid_t *grid;
2591     int                 i, j4, j, si, b;
2592     int                 c[GPU_NSUBCELL+1];
2593
2594     /* This code only produces correct statistics with domain decomposition */
2595     grid = &nbs->grid[0];
2596
2597     fprintf(fp, "nbl nsci %d ncj4 %d nsi %d excl4 %d\n",
2598             nbl->nsci, nbl->ncj4, nbl->nci_tot, nbl->nexcl);
2599     fprintf(fp, "nbl na_c %d rl %g ncp %d per cell %.1f atoms %.1f ratio %.2f\n",
2600             nbl->na_ci, rl, nbl->nci_tot, nbl->nci_tot/(double)grid->nsubc_tot,
2601             nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c,
2602             nbl->nci_tot/(double)grid->nsubc_tot*grid->na_c/(0.5*4.0/3.0*M_PI*rl*rl*rl*grid->nsubc_tot*grid->na_c/det(nbs->box)));
2603
2604     fprintf(fp, "nbl average j super cell list length %.1f\n",
2605             0.25*nbl->ncj4/(double)nbl->nsci);
2606     fprintf(fp, "nbl average i sub cell list length %.1f\n",
2607             nbl->nci_tot/((double)nbl->ncj4));
2608
2609     for (si = 0; si <= GPU_NSUBCELL; si++)
2610     {
2611         c[si] = 0;
2612     }
2613     for (i = 0; i < nbl->nsci; i++)
2614     {
2615         for (j4 = nbl->sci[i].cj4_ind_start; j4 < nbl->sci[i].cj4_ind_end; j4++)
2616         {
2617             for (j = 0; j < NBNXN_GPU_JGROUP_SIZE; j++)
2618             {
2619                 b = 0;
2620                 for (si = 0; si < GPU_NSUBCELL; si++)
2621                 {
2622                     if (nbl->cj4[j4].imei[0].imask & (1U << (j*GPU_NSUBCELL + si)))
2623                     {
2624                         b++;
2625                     }
2626                 }
2627                 c[b]++;
2628             }
2629         }
2630     }
2631     for (b = 0; b <= GPU_NSUBCELL; b++)
2632     {
2633         fprintf(fp, "nbl j-list #i-subcell %d %7d %4.1f\n",
2634                 b, c[b], 100.0*c[b]/(double)(nbl->ncj4*NBNXN_GPU_JGROUP_SIZE));
2635     }
2636 }
2637
2638 /* Returns a pointer to the exclusion mask for cj4-unit cj4, warp warp */
2639 static void low_get_nbl_exclusions(nbnxn_pairlist_t *nbl, int cj4,
2640                                    int warp, nbnxn_excl_t **excl)
2641 {
2642     if (nbl->cj4[cj4].imei[warp].excl_ind == 0)
2643     {
2644         /* No exclusions set, make a new list entry */
2645         nbl->cj4[cj4].imei[warp].excl_ind = nbl->nexcl;
2646         nbl->nexcl++;
2647         *excl = &nbl->excl[nbl->cj4[cj4].imei[warp].excl_ind];
2648         set_no_excls(*excl);
2649     }
2650     else
2651     {
2652         /* We already have some exclusions, new ones can be added to the list */
2653         *excl = &nbl->excl[nbl->cj4[cj4].imei[warp].excl_ind];
2654     }
2655 }
2656
2657 /* Returns a pointer to the exclusion mask for cj4-unit cj4, warp warp,
2658  * allocates extra memory, if necessary.
2659  */
2660 static void get_nbl_exclusions_1(nbnxn_pairlist_t *nbl, int cj4,
2661                                  int warp, nbnxn_excl_t **excl)
2662 {
2663     if (nbl->cj4[cj4].imei[warp].excl_ind == 0)
2664     {
2665         /* We need to make a new list entry, check if we have space */
2666         check_excl_space(nbl, 1);
2667     }
2668     low_get_nbl_exclusions(nbl, cj4, warp, excl);
2669 }
2670
2671 /* Returns pointers to the exclusion mask for cj4-unit cj4 for both warps,
2672  * allocates extra memory, if necessary.
2673  */
2674 static void get_nbl_exclusions_2(nbnxn_pairlist_t *nbl, int cj4,
2675                                  nbnxn_excl_t **excl_w0,
2676                                  nbnxn_excl_t **excl_w1)
2677 {
2678     /* Check for space we might need */
2679     check_excl_space(nbl, 2);
2680
2681     low_get_nbl_exclusions(nbl, cj4, 0, excl_w0);
2682     low_get_nbl_exclusions(nbl, cj4, 1, excl_w1);
2683 }
2684
2685 /* Sets the self exclusions i=j and pair exclusions i>j */
2686 static void set_self_and_newton_excls_supersub(nbnxn_pairlist_t *nbl,
2687                                                int cj4_ind, int sj_offset,
2688                                                int si)
2689 {
2690     nbnxn_excl_t *excl[2];
2691     int           ei, ej, w;
2692
2693     /* Here we only set the set self and double pair exclusions */
2694
2695     get_nbl_exclusions_2(nbl, cj4_ind, &excl[0], &excl[1]);
2696
2697     /* Only minor < major bits set */
2698     for (ej = 0; ej < nbl->na_ci; ej++)
2699     {
2700         w = (ej>>2);
2701         for (ei = ej; ei < nbl->na_ci; ei++)
2702         {
2703             excl[w]->pair[(ej & (NBNXN_GPU_JGROUP_SIZE-1))*nbl->na_ci + ei] &=
2704                 ~(1U << (sj_offset*GPU_NSUBCELL + si));
2705         }
2706     }
2707 }
2708
2709 /* Returns a diagonal or off-diagonal interaction mask for plain C lists */
2710 static unsigned int get_imask(gmx_bool rdiag, int ci, int cj)
2711 {
2712     return (rdiag && ci == cj ? NBNXN_INTERACTION_MASK_DIAG : NBNXN_INTERACTION_MASK_ALL);
2713 }
2714
2715 /* Returns a diagonal or off-diagonal interaction mask for cj-size=2 */
2716 static unsigned int get_imask_simd_j2(gmx_bool rdiag, int ci, int cj)
2717 {
2718     return (rdiag && ci*2 == cj ? NBNXN_INTERACTION_MASK_DIAG_J2_0 :
2719             (rdiag && ci*2+1 == cj ? NBNXN_INTERACTION_MASK_DIAG_J2_1 :
2720              NBNXN_INTERACTION_MASK_ALL));
2721 }
2722
2723 /* Returns a diagonal or off-diagonal interaction mask for cj-size=4 */
2724 static unsigned int get_imask_simd_j4(gmx_bool rdiag, int ci, int cj)
2725 {
2726     return (rdiag && ci == cj ? NBNXN_INTERACTION_MASK_DIAG : NBNXN_INTERACTION_MASK_ALL);
2727 }
2728
2729 /* Returns a diagonal or off-diagonal interaction mask for cj-size=8 */
2730 static unsigned int get_imask_simd_j8(gmx_bool rdiag, int ci, int cj)
2731 {
2732     return (rdiag && ci == cj*2 ? NBNXN_INTERACTION_MASK_DIAG_J8_0 :
2733             (rdiag && ci == cj*2+1 ? NBNXN_INTERACTION_MASK_DIAG_J8_1 :
2734              NBNXN_INTERACTION_MASK_ALL));
2735 }
2736
2737 #ifdef GMX_NBNXN_SIMD
2738 #if GMX_SIMD_WIDTH_HERE == 2
2739 #define get_imask_simd_4xn  get_imask_simd_j2
2740 #endif
2741 #if GMX_SIMD_WIDTH_HERE == 4
2742 #define get_imask_simd_4xn  get_imask_simd_j4
2743 #endif
2744 #if GMX_SIMD_WIDTH_HERE == 8
2745 #define get_imask_simd_4xn  get_imask_simd_j8
2746 #define get_imask_simd_2xnn get_imask_simd_j4
2747 #endif
2748 #if GMX_SIMD_WIDTH_HERE == 16
2749 #define get_imask_simd_2xnn get_imask_simd_j8
2750 #endif
2751 #endif
2752
2753 /* Plain C code for making a pair list of cell ci vs cell cjf-cjl.
2754  * Checks bounding box distances and possibly atom pair distances.
2755  */
2756 static void make_cluster_list_simple(const nbnxn_grid_t *gridj,
2757                                      nbnxn_pairlist_t *nbl,
2758                                      int ci, int cjf, int cjl,
2759                                      gmx_bool remove_sub_diag,
2760                                      const real *x_j,
2761                                      real rl2, float rbb2,
2762                                      int *ndistc)
2763 {
2764     const nbnxn_list_work_t *work;
2765
2766     const nbnxn_bb_t        *bb_ci;
2767     const real              *x_ci;
2768
2769     gmx_bool                 InRange;
2770     real                     d2;
2771     int                      cjf_gl, cjl_gl, cj;
2772
2773     work = nbl->work;
2774
2775     bb_ci = nbl->work->bb_ci;
2776     x_ci  = nbl->work->x_ci;
2777
2778     InRange = FALSE;
2779     while (!InRange && cjf <= cjl)
2780     {
2781         d2       = subc_bb_dist2(0, bb_ci, cjf, gridj->bb);
2782         *ndistc += 2;
2783
2784         /* Check if the distance is within the distance where
2785          * we use only the bounding box distance rbb,
2786          * or within the cut-off and there is at least one atom pair
2787          * within the cut-off.
2788          */
2789         if (d2 < rbb2)
2790         {
2791             InRange = TRUE;
2792         }
2793         else if (d2 < rl2)
2794         {
2795             int i, j;
2796
2797             cjf_gl = gridj->cell0 + cjf;
2798             for (i = 0; i < NBNXN_CPU_CLUSTER_I_SIZE && !InRange; i++)
2799             {
2800                 for (j = 0; j < NBNXN_CPU_CLUSTER_I_SIZE; j++)
2801                 {
2802                     InRange = InRange ||
2803                         (sqr(x_ci[i*STRIDE_XYZ+XX] - x_j[(cjf_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+XX]) +
2804                          sqr(x_ci[i*STRIDE_XYZ+YY] - x_j[(cjf_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+YY]) +
2805                          sqr(x_ci[i*STRIDE_XYZ+ZZ] - x_j[(cjf_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+ZZ]) < rl2);
2806                 }
2807             }
2808             *ndistc += NBNXN_CPU_CLUSTER_I_SIZE*NBNXN_CPU_CLUSTER_I_SIZE;
2809         }
2810         if (!InRange)
2811         {
2812             cjf++;
2813         }
2814     }
2815     if (!InRange)
2816     {
2817         return;
2818     }
2819
2820     InRange = FALSE;
2821     while (!InRange && cjl > cjf)
2822     {
2823         d2       = subc_bb_dist2(0, bb_ci, cjl, gridj->bb);
2824         *ndistc += 2;
2825
2826         /* Check if the distance is within the distance where
2827          * we use only the bounding box distance rbb,
2828          * or within the cut-off and there is at least one atom pair
2829          * within the cut-off.
2830          */
2831         if (d2 < rbb2)
2832         {
2833             InRange = TRUE;
2834         }
2835         else if (d2 < rl2)
2836         {
2837             int i, j;
2838
2839             cjl_gl = gridj->cell0 + cjl;
2840             for (i = 0; i < NBNXN_CPU_CLUSTER_I_SIZE && !InRange; i++)
2841             {
2842                 for (j = 0; j < NBNXN_CPU_CLUSTER_I_SIZE; j++)
2843                 {
2844                     InRange = InRange ||
2845                         (sqr(x_ci[i*STRIDE_XYZ+XX] - x_j[(cjl_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+XX]) +
2846                          sqr(x_ci[i*STRIDE_XYZ+YY] - x_j[(cjl_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+YY]) +
2847                          sqr(x_ci[i*STRIDE_XYZ+ZZ] - x_j[(cjl_gl*NBNXN_CPU_CLUSTER_I_SIZE+j)*STRIDE_XYZ+ZZ]) < rl2);
2848                 }
2849             }
2850             *ndistc += NBNXN_CPU_CLUSTER_I_SIZE*NBNXN_CPU_CLUSTER_I_SIZE;
2851         }
2852         if (!InRange)
2853         {
2854             cjl--;
2855         }
2856     }
2857
2858     if (cjf <= cjl)
2859     {
2860         for (cj = cjf; cj <= cjl; cj++)
2861         {
2862             /* Store cj and the interaction mask */
2863             nbl->cj[nbl->ncj].cj   = gridj->cell0 + cj;
2864             nbl->cj[nbl->ncj].excl = get_imask(remove_sub_diag, ci, cj);
2865             nbl->ncj++;
2866         }
2867         /* Increase the closing index in i super-cell list */
2868         nbl->ci[nbl->nci].cj_ind_end = nbl->ncj;
2869     }
2870 }
2871
2872 #ifdef GMX_NBNXN_SIMD_4XN
2873 #include "nbnxn_search_simd_4xn.h"
2874 #endif
2875 #ifdef GMX_NBNXN_SIMD_2XNN
2876 #include "nbnxn_search_simd_2xnn.h"
2877 #endif
2878
2879 /* Plain C or SIMD4 code for making a pair list of super-cell sci vs scj.
2880  * Checks bounding box distances and possibly atom pair distances.
2881  */
2882 static void make_cluster_list_supersub(const nbnxn_grid_t *gridi,
2883                                        const nbnxn_grid_t *gridj,
2884                                        nbnxn_pairlist_t *nbl,
2885                                        int sci, int scj,
2886                                        gmx_bool sci_equals_scj,
2887                                        int stride, const real *x,
2888                                        real rl2, float rbb2,
2889                                        int *ndistc)
2890 {
2891     int               na_c;
2892     int               npair;
2893     int               cjo, ci1, ci, cj, cj_gl;
2894     int               cj4_ind, cj_offset;
2895     unsigned          imask;
2896     nbnxn_cj4_t      *cj4;
2897 #ifdef NBNXN_BBXXXX
2898     const float      *pbb_ci;
2899 #else
2900     const nbnxn_bb_t *bb_ci;
2901 #endif
2902     const real       *x_ci;
2903     float            *d2l, d2;
2904     int               w;
2905 #define PRUNE_LIST_CPU_ONE
2906 #ifdef PRUNE_LIST_CPU_ONE
2907     int  ci_last = -1;
2908 #endif
2909
2910     d2l = nbl->work->d2;
2911
2912 #ifdef NBNXN_BBXXXX
2913     pbb_ci = nbl->work->pbb_ci;
2914 #else
2915     bb_ci  = nbl->work->bb_ci;
2916 #endif
2917     x_ci   = nbl->work->x_ci;
2918
2919     na_c = gridj->na_c;
2920
2921     for (cjo = 0; cjo < gridj->nsubc[scj]; cjo++)
2922     {
2923         cj4_ind   = (nbl->work->cj_ind >> NBNXN_GPU_JGROUP_SIZE_2LOG);
2924         cj_offset = nbl->work->cj_ind - cj4_ind*NBNXN_GPU_JGROUP_SIZE;
2925         cj4       = &nbl->cj4[cj4_ind];
2926
2927         cj = scj*GPU_NSUBCELL + cjo;
2928
2929         cj_gl = gridj->cell0*GPU_NSUBCELL + cj;
2930
2931         /* Initialize this j-subcell i-subcell list */
2932         cj4->cj[cj_offset] = cj_gl;
2933         imask              = 0;
2934
2935         if (sci_equals_scj)
2936         {
2937             ci1 = cjo + 1;
2938         }
2939         else
2940         {
2941             ci1 = gridi->nsubc[sci];
2942         }
2943
2944 #ifdef NBNXN_BBXXXX
2945         /* Determine all ci1 bb distances in one call with SIMD4 */
2946         subc_bb_dist2_simd4_xxxx(gridj->pbb+(cj>>STRIDE_PBB_2LOG)*NNBSBB_XXXX+(cj & (STRIDE_PBB-1)),
2947                                  ci1, pbb_ci, d2l);
2948         *ndistc += na_c*2;
2949 #endif
2950
2951         npair = 0;
2952         /* We use a fixed upper-bound instead of ci1 to help optimization */
2953         for (ci = 0; ci < GPU_NSUBCELL; ci++)
2954         {
2955             if (ci == ci1)
2956             {
2957                 break;
2958             }
2959
2960 #ifndef NBNXN_BBXXXX
2961             /* Determine the bb distance between ci and cj */
2962             d2l[ci]  = subc_bb_dist2(ci, bb_ci, cj, gridj->bb);
2963             *ndistc += 2;
2964 #endif
2965             d2 = d2l[ci];
2966
2967 #ifdef PRUNE_LIST_CPU_ALL
2968             /* Check if the distance is within the distance where
2969              * we use only the bounding box distance rbb,
2970              * or within the cut-off and there is at least one atom pair
2971              * within the cut-off. This check is very costly.
2972              */
2973             *ndistc += na_c*na_c;
2974             if (d2 < rbb2 ||
2975                 (d2 < rl2 &&
2976 #ifdef NBNXN_PBB_SIMD4
2977                  subc_in_range_simd4
2978 #else
2979                  subc_in_range_x
2980 #endif
2981                      (na_c, ci, x_ci, cj_gl, stride, x, rl2)))
2982 #else
2983             /* Check if the distance between the two bounding boxes
2984              * in within the pair-list cut-off.
2985              */
2986             if (d2 < rl2)
2987 #endif
2988             {
2989                 /* Flag this i-subcell to be taken into account */
2990                 imask |= (1U << (cj_offset*GPU_NSUBCELL+ci));
2991
2992 #ifdef PRUNE_LIST_CPU_ONE
2993                 ci_last = ci;
2994 #endif
2995
2996                 npair++;
2997             }
2998         }
2999
3000 #ifdef PRUNE_LIST_CPU_ONE
3001         /* If we only found 1 pair, check if any atoms are actually
3002          * within the cut-off, so we could get rid of it.
3003          */
3004         if (npair == 1 && d2l[ci_last] >= rbb2)
3005         {
3006             /* Avoid using function pointers here, as it's slower */
3007             if (
3008 #ifdef NBNXN_PBB_SIMD4
3009                 !subc_in_range_simd4
3010 #else
3011                 !subc_in_range_x
3012 #endif
3013                     (na_c, ci_last, x_ci, cj_gl, stride, x, rl2))
3014             {
3015                 imask &= ~(1U << (cj_offset*GPU_NSUBCELL+ci_last));
3016                 npair--;
3017             }
3018         }
3019 #endif
3020
3021         if (npair > 0)
3022         {
3023             /* We have a useful sj entry, close it now */
3024
3025             /* Set the exclucions for the ci== sj entry.
3026              * Here we don't bother to check if this entry is actually flagged,
3027              * as it will nearly always be in the list.
3028              */
3029             if (sci_equals_scj)
3030             {
3031                 set_self_and_newton_excls_supersub(nbl, cj4_ind, cj_offset, cjo);
3032             }
3033
3034             /* Copy the cluster interaction mask to the list */
3035             for (w = 0; w < NWARP; w++)
3036             {
3037                 cj4->imei[w].imask |= imask;
3038             }
3039
3040             nbl->work->cj_ind++;
3041
3042             /* Keep the count */
3043             nbl->nci_tot += npair;
3044
3045             /* Increase the closing index in i super-cell list */
3046             nbl->sci[nbl->nsci].cj4_ind_end =
3047                 ((nbl->work->cj_ind+NBNXN_GPU_JGROUP_SIZE-1) >> NBNXN_GPU_JGROUP_SIZE_2LOG);
3048         }
3049     }
3050 }
3051
3052 /* Set all atom-pair exclusions from the topology stored in excl
3053  * as masks in the pair-list for simple list i-entry nbl_ci
3054  */
3055 static void set_ci_top_excls(const nbnxn_search_t nbs,
3056                              nbnxn_pairlist_t    *nbl,
3057                              gmx_bool             diagRemoved,
3058                              int                  na_ci_2log,
3059                              int                  na_cj_2log,
3060                              const nbnxn_ci_t    *nbl_ci,
3061                              const t_blocka      *excl)
3062 {
3063     const int    *cell;
3064     int           ci;
3065     int           cj_ind_first, cj_ind_last;
3066     int           cj_first, cj_last;
3067     int           ndirect;
3068     int           i, ai, aj, si, eind, ge, se;
3069     int           found, cj_ind_0, cj_ind_1, cj_ind_m;
3070     int           cj_m;
3071     gmx_bool      Found_si;
3072     int           si_ind;
3073     nbnxn_excl_t *nbl_excl;
3074     int           inner_i, inner_e;
3075
3076     cell = nbs->cell;
3077
3078     if (nbl_ci->cj_ind_end == nbl_ci->cj_ind_start)
3079     {
3080         /* Empty list */
3081         return;
3082     }
3083
3084     ci = nbl_ci->ci;
3085
3086     cj_ind_first = nbl_ci->cj_ind_start;
3087     cj_ind_last  = nbl->ncj - 1;
3088
3089     cj_first = nbl->cj[cj_ind_first].cj;
3090     cj_last  = nbl->cj[cj_ind_last].cj;
3091
3092     /* Determine how many contiguous j-cells we have starting
3093      * from the first i-cell. This number can be used to directly
3094      * calculate j-cell indices for excluded atoms.
3095      */
3096     ndirect = 0;
3097     if (na_ci_2log == na_cj_2log)
3098     {
3099         while (cj_ind_first + ndirect <= cj_ind_last &&
3100                nbl->cj[cj_ind_first+ndirect].cj == ci + ndirect)
3101         {
3102             ndirect++;
3103         }
3104     }
3105 #ifdef NBNXN_SEARCH_BB_SIMD4
3106     else
3107     {
3108         while (cj_ind_first + ndirect <= cj_ind_last &&
3109                nbl->cj[cj_ind_first+ndirect].cj == ci_to_cj(na_cj_2log, ci) + ndirect)
3110         {
3111             ndirect++;
3112         }
3113     }
3114 #endif
3115
3116     /* Loop over the atoms in the i super-cell */
3117     for (i = 0; i < nbl->na_sc; i++)
3118     {
3119         ai = nbs->a[ci*nbl->na_sc+i];
3120         if (ai >= 0)
3121         {
3122             si  = (i>>na_ci_2log);
3123
3124             /* Loop over the topology-based exclusions for this i-atom */
3125             for (eind = excl->index[ai]; eind < excl->index[ai+1]; eind++)
3126             {
3127                 aj = excl->a[eind];
3128
3129                 if (aj == ai)
3130                 {
3131                     /* The self exclusion are already set, save some time */
3132                     continue;
3133                 }
3134
3135                 ge = cell[aj];
3136
3137                 /* Without shifts we only calculate interactions j>i
3138                  * for one-way pair-lists.
3139                  */
3140                 if (diagRemoved && ge <= ci*nbl->na_sc + i)
3141                 {
3142                     continue;
3143                 }
3144
3145                 se = (ge >> na_cj_2log);
3146
3147                 /* Could the cluster se be in our list? */
3148                 if (se >= cj_first && se <= cj_last)
3149                 {
3150                     if (se < cj_first + ndirect)
3151                     {
3152                         /* We can calculate cj_ind directly from se */
3153                         found = cj_ind_first + se - cj_first;
3154                     }
3155                     else
3156                     {
3157                         /* Search for se using bisection */
3158                         found    = -1;
3159                         cj_ind_0 = cj_ind_first + ndirect;
3160                         cj_ind_1 = cj_ind_last + 1;
3161                         while (found == -1 && cj_ind_0 < cj_ind_1)
3162                         {
3163                             cj_ind_m = (cj_ind_0 + cj_ind_1)>>1;
3164
3165                             cj_m = nbl->cj[cj_ind_m].cj;
3166
3167                             if (se == cj_m)
3168                             {
3169                                 found = cj_ind_m;
3170                             }
3171                             else if (se < cj_m)
3172                             {
3173                                 cj_ind_1 = cj_ind_m;
3174                             }
3175                             else
3176                             {
3177                                 cj_ind_0 = cj_ind_m + 1;
3178                             }
3179                         }
3180                     }
3181
3182                     if (found >= 0)
3183                     {
3184                         inner_i = i  - (si << na_ci_2log);
3185                         inner_e = ge - (se << na_cj_2log);
3186
3187                         nbl->cj[found].excl &= ~(1U<<((inner_i<<na_cj_2log) + inner_e));
3188 /* The next code line is usually not needed. We do not want to version
3189  * away the above line, because there is logic that relies on being
3190  * able to detect easily whether any exclusions exist. */
3191 #if (defined GMX_CPU_ACCELERATION_IBM_QPX)
3192                         nbl->cj[found].interaction_mask_indices[inner_i] &= ~(1U << inner_e);
3193 #endif
3194                     }
3195                 }
3196             }
3197         }
3198     }
3199 }
3200
3201 /* Set all atom-pair exclusions from the topology stored in excl
3202  * as masks in the pair-list for i-super-cell entry nbl_sci
3203  */
3204 static void set_sci_top_excls(const nbnxn_search_t nbs,
3205                               nbnxn_pairlist_t    *nbl,
3206                               gmx_bool             diagRemoved,
3207                               int                  na_c_2log,
3208                               const nbnxn_sci_t   *nbl_sci,
3209                               const t_blocka      *excl)
3210 {
3211     const int    *cell;
3212     int           na_c;
3213     int           sci;
3214     int           cj_ind_first, cj_ind_last;
3215     int           cj_first, cj_last;
3216     int           ndirect;
3217     int           i, ai, aj, si, eind, ge, se;
3218     int           found, cj_ind_0, cj_ind_1, cj_ind_m;
3219     int           cj_m;
3220     gmx_bool      Found_si;
3221     int           si_ind;
3222     nbnxn_excl_t *nbl_excl;
3223     int           inner_i, inner_e, w;
3224
3225     cell = nbs->cell;
3226
3227     na_c = nbl->na_ci;
3228
3229     if (nbl_sci->cj4_ind_end == nbl_sci->cj4_ind_start)
3230     {
3231         /* Empty list */
3232         return;
3233     }
3234
3235     sci = nbl_sci->sci;
3236
3237     cj_ind_first = nbl_sci->cj4_ind_start*NBNXN_GPU_JGROUP_SIZE;
3238     cj_ind_last  = nbl->work->cj_ind - 1;
3239
3240     cj_first = nbl->cj4[nbl_sci->cj4_ind_start].cj[0];
3241     cj_last  = nbl_cj(nbl, cj_ind_last);
3242
3243     /* Determine how many contiguous j-clusters we have starting
3244      * from the first i-cluster. This number can be used to directly
3245      * calculate j-cluster indices for excluded atoms.
3246      */
3247     ndirect = 0;
3248     while (cj_ind_first + ndirect <= cj_ind_last &&
3249            nbl_cj(nbl, cj_ind_first+ndirect) == sci*GPU_NSUBCELL + ndirect)
3250     {
3251         ndirect++;
3252     }
3253
3254     /* Loop over the atoms in the i super-cell */
3255     for (i = 0; i < nbl->na_sc; i++)
3256     {
3257         ai = nbs->a[sci*nbl->na_sc+i];
3258         if (ai >= 0)
3259         {
3260             si  = (i>>na_c_2log);
3261
3262             /* Loop over the topology-based exclusions for this i-atom */
3263             for (eind = excl->index[ai]; eind < excl->index[ai+1]; eind++)
3264             {
3265                 aj = excl->a[eind];
3266
3267                 if (aj == ai)
3268                 {
3269                     /* The self exclusion are already set, save some time */
3270                     continue;
3271                 }
3272
3273                 ge = cell[aj];
3274
3275                 /* Without shifts we only calculate interactions j>i
3276                  * for one-way pair-lists.
3277                  */
3278                 if (diagRemoved && ge <= sci*nbl->na_sc + i)
3279                 {
3280                     continue;
3281                 }
3282
3283                 se = ge>>na_c_2log;
3284                 /* Could the cluster se be in our list? */
3285                 if (se >= cj_first && se <= cj_last)
3286                 {
3287                     if (se < cj_first + ndirect)
3288                     {
3289                         /* We can calculate cj_ind directly from se */
3290                         found = cj_ind_first + se - cj_first;
3291                     }
3292                     else
3293                     {
3294                         /* Search for se using bisection */
3295                         found    = -1;
3296                         cj_ind_0 = cj_ind_first + ndirect;
3297                         cj_ind_1 = cj_ind_last + 1;
3298                         while (found == -1 && cj_ind_0 < cj_ind_1)
3299                         {
3300                             cj_ind_m = (cj_ind_0 + cj_ind_1)>>1;
3301
3302                             cj_m = nbl_cj(nbl, cj_ind_m);
3303
3304                             if (se == cj_m)
3305                             {
3306                                 found = cj_ind_m;
3307                             }
3308                             else if (se < cj_m)
3309                             {
3310                                 cj_ind_1 = cj_ind_m;
3311                             }
3312                             else
3313                             {
3314                                 cj_ind_0 = cj_ind_m + 1;
3315                             }
3316                         }
3317                     }
3318
3319                     if (found >= 0)
3320                     {
3321                         inner_i = i  - si*na_c;
3322                         inner_e = ge - se*na_c;
3323
3324 /* Macro for getting the index of atom a within a cluster */
3325 #define AMODCJ4(a)  ((a) & (NBNXN_GPU_JGROUP_SIZE - 1))
3326 /* Macro for converting an atom number to a cluster number */
3327 #define A2CJ4(a)    ((a) >> NBNXN_GPU_JGROUP_SIZE_2LOG)
3328 /* Macro for getting the index of an i-atom within a warp */
3329 #define AMODWI(a)   ((a) & (NBNXN_GPU_CLUSTER_SIZE/2 - 1))
3330
3331                         if (nbl_imask0(nbl, found) & (1U << (AMODCJ4(found)*GPU_NSUBCELL + si)))
3332                         {
3333                             w       = (inner_e >> 2);
3334
3335                             get_nbl_exclusions_1(nbl, A2CJ4(found), w, &nbl_excl);
3336
3337                             nbl_excl->pair[AMODWI(inner_e)*nbl->na_ci+inner_i] &=
3338                                 ~(1U << (AMODCJ4(found)*GPU_NSUBCELL + si));
3339                         }
3340
3341 #undef AMODCJ4
3342 #undef A2CJ4
3343 #undef AMODWI
3344                     }
3345                 }
3346             }
3347         }
3348     }
3349 }
3350
3351 /* Reallocate the simple ci list for at least n entries */
3352 static void nb_realloc_ci(nbnxn_pairlist_t *nbl, int n)
3353 {
3354     nbl->ci_nalloc = over_alloc_small(n);
3355     nbnxn_realloc_void((void **)&nbl->ci,
3356                        nbl->nci*sizeof(*nbl->ci),
3357                        nbl->ci_nalloc*sizeof(*nbl->ci),
3358                        nbl->alloc, nbl->free);
3359 }
3360
3361 /* Reallocate the super-cell sci list for at least n entries */
3362 static void nb_realloc_sci(nbnxn_pairlist_t *nbl, int n)
3363 {
3364     nbl->sci_nalloc = over_alloc_small(n);
3365     nbnxn_realloc_void((void **)&nbl->sci,
3366                        nbl->nsci*sizeof(*nbl->sci),
3367                        nbl->sci_nalloc*sizeof(*nbl->sci),
3368                        nbl->alloc, nbl->free);
3369 }
3370
3371 /* Make a new ci entry at index nbl->nci */
3372 static void new_ci_entry(nbnxn_pairlist_t *nbl, int ci, int shift, int flags)
3373 {
3374     if (nbl->nci + 1 > nbl->ci_nalloc)
3375     {
3376         nb_realloc_ci(nbl, nbl->nci+1);
3377     }
3378     nbl->ci[nbl->nci].ci            = ci;
3379     nbl->ci[nbl->nci].shift         = shift;
3380     /* Store the interaction flags along with the shift */
3381     nbl->ci[nbl->nci].shift        |= flags;
3382     nbl->ci[nbl->nci].cj_ind_start  = nbl->ncj;
3383     nbl->ci[nbl->nci].cj_ind_end    = nbl->ncj;
3384 }
3385
3386 /* Make a new sci entry at index nbl->nsci */
3387 static void new_sci_entry(nbnxn_pairlist_t *nbl, int sci, int shift)
3388 {
3389     if (nbl->nsci + 1 > nbl->sci_nalloc)
3390     {
3391         nb_realloc_sci(nbl, nbl->nsci+1);
3392     }
3393     nbl->sci[nbl->nsci].sci           = sci;
3394     nbl->sci[nbl->nsci].shift         = shift;
3395     nbl->sci[nbl->nsci].cj4_ind_start = nbl->ncj4;
3396     nbl->sci[nbl->nsci].cj4_ind_end   = nbl->ncj4;
3397 }
3398
3399 /* Sort the simple j-list cj on exclusions.
3400  * Entries with exclusions will all be sorted to the beginning of the list.
3401  */
3402 static void sort_cj_excl(nbnxn_cj_t *cj, int ncj,
3403                          nbnxn_list_work_t *work)
3404 {
3405     int jnew, j;
3406
3407     if (ncj > work->cj_nalloc)
3408     {
3409         work->cj_nalloc = over_alloc_large(ncj);
3410         srenew(work->cj, work->cj_nalloc);
3411     }
3412
3413     /* Make a list of the j-cells involving exclusions */
3414     jnew = 0;
3415     for (j = 0; j < ncj; j++)
3416     {
3417         if (cj[j].excl != NBNXN_INTERACTION_MASK_ALL)
3418         {
3419             work->cj[jnew++] = cj[j];
3420         }
3421     }
3422     /* Check if there are exclusions at all or not just the first entry */
3423     if (!((jnew == 0) ||
3424           (jnew == 1 && cj[0].excl != NBNXN_INTERACTION_MASK_ALL)))
3425     {
3426         for (j = 0; j < ncj; j++)
3427         {
3428             if (cj[j].excl == NBNXN_INTERACTION_MASK_ALL)
3429             {
3430                 work->cj[jnew++] = cj[j];
3431             }
3432         }
3433         for (j = 0; j < ncj; j++)
3434         {
3435             cj[j] = work->cj[j];
3436         }
3437     }
3438 }
3439
3440 /* Close this simple list i entry */
3441 static void close_ci_entry_simple(nbnxn_pairlist_t *nbl)
3442 {
3443     int jlen;
3444
3445     /* All content of the new ci entry have already been filled correctly,
3446      * we only need to increase the count here (for non empty lists).
3447      */
3448     jlen = nbl->ci[nbl->nci].cj_ind_end - nbl->ci[nbl->nci].cj_ind_start;
3449     if (jlen > 0)
3450     {
3451         sort_cj_excl(nbl->cj+nbl->ci[nbl->nci].cj_ind_start, jlen, nbl->work);
3452
3453         /* The counts below are used for non-bonded pair/flop counts
3454          * and should therefore match the available kernel setups.
3455          */
3456         if (!(nbl->ci[nbl->nci].shift & NBNXN_CI_DO_COUL(0)))
3457         {
3458             nbl->work->ncj_noq += jlen;
3459         }
3460         else if ((nbl->ci[nbl->nci].shift & NBNXN_CI_HALF_LJ(0)) ||
3461                  !(nbl->ci[nbl->nci].shift & NBNXN_CI_DO_LJ(0)))
3462         {
3463             nbl->work->ncj_hlj += jlen;
3464         }
3465
3466         nbl->nci++;
3467     }
3468 }
3469
3470 /* Split sci entry for load balancing on the GPU.
3471  * Splitting ensures we have enough lists to fully utilize the whole GPU.
3472  * With progBal we generate progressively smaller lists, which improves
3473  * load balancing. As we only know the current count on our own thread,
3474  * we will need to estimate the current total amount of i-entries.
3475  * As the lists get concatenated later, this estimate depends
3476  * both on nthread and our own thread index.
3477  */
3478 static void split_sci_entry(nbnxn_pairlist_t *nbl,
3479                             int nsp_max_av, gmx_bool progBal, int nc_bal,
3480                             int thread, int nthread)
3481 {
3482     int nsci_est;
3483     int nsp_max;
3484     int cj4_start, cj4_end, j4len, cj4;
3485     int sci;
3486     int nsp, nsp_sci, nsp_cj4, nsp_cj4_e, nsp_cj4_p;
3487     int p;
3488
3489     if (progBal)
3490     {
3491         /* Estimate the total numbers of ci's of the nblist combined
3492          * over all threads using the target number of ci's.
3493          */
3494         nsci_est = nc_bal*thread/nthread + nbl->nsci;
3495
3496         /* The first ci blocks should be larger, to avoid overhead.
3497          * The last ci blocks should be smaller, to improve load balancing.
3498          */
3499         nsp_max = max(1,
3500                       nsp_max_av*nc_bal*3/(2*(nsci_est - 1 + nc_bal)));
3501     }
3502     else
3503     {
3504         nsp_max = nsp_max_av;
3505     }
3506
3507     cj4_start = nbl->sci[nbl->nsci-1].cj4_ind_start;
3508     cj4_end   = nbl->sci[nbl->nsci-1].cj4_ind_end;
3509     j4len     = cj4_end - cj4_start;
3510
3511     if (j4len > 1 && j4len*GPU_NSUBCELL*NBNXN_GPU_JGROUP_SIZE > nsp_max)
3512     {
3513         /* Remove the last ci entry and process the cj4's again */
3514         nbl->nsci -= 1;
3515
3516         sci        = nbl->nsci;
3517         nsp        = 0;
3518         nsp_sci    = 0;
3519         nsp_cj4_e  = 0;
3520         nsp_cj4    = 0;
3521         for (cj4 = cj4_start; cj4 < cj4_end; cj4++)
3522         {
3523             nsp_cj4_p = nsp_cj4;
3524             /* Count the number of cluster pairs in this cj4 group */
3525             nsp_cj4   = 0;
3526             for (p = 0; p < GPU_NSUBCELL*NBNXN_GPU_JGROUP_SIZE; p++)
3527             {
3528                 nsp_cj4 += (nbl->cj4[cj4].imei[0].imask >> p) & 1;
3529             }
3530
3531             if (nsp_cj4 > 0 && nsp + nsp_cj4 > nsp_max)
3532             {
3533                 /* Split the list at cj4 */
3534                 nbl->sci[sci].cj4_ind_end = cj4;
3535                 /* Create a new sci entry */
3536                 sci++;
3537                 nbl->nsci++;
3538                 if (nbl->nsci+1 > nbl->sci_nalloc)
3539                 {
3540                     nb_realloc_sci(nbl, nbl->nsci+1);
3541                 }
3542                 nbl->sci[sci].sci           = nbl->sci[nbl->nsci-1].sci;
3543                 nbl->sci[sci].shift         = nbl->sci[nbl->nsci-1].shift;
3544                 nbl->sci[sci].cj4_ind_start = cj4;
3545                 nsp_sci                     = nsp;
3546                 nsp_cj4_e                   = nsp_cj4_p;
3547                 nsp                         = 0;
3548             }
3549             nsp += nsp_cj4;
3550         }
3551
3552         /* Put the remaining cj4's in the last sci entry */
3553         nbl->sci[sci].cj4_ind_end = cj4_end;
3554
3555         /* Possibly balance out the last two sci's
3556          * by moving the last cj4 of the second last sci.
3557          */
3558         if (nsp_sci - nsp_cj4_e >= nsp + nsp_cj4_e)
3559         {
3560             nbl->sci[sci-1].cj4_ind_end--;
3561             nbl->sci[sci].cj4_ind_start--;
3562         }
3563
3564         nbl->nsci++;
3565     }
3566 }
3567
3568 /* Clost this super/sub list i entry */
3569 static void close_ci_entry_supersub(nbnxn_pairlist_t *nbl,
3570                                     int nsp_max_av,
3571                                     gmx_bool progBal, int nc_bal,
3572                                     int thread, int nthread)
3573 {
3574     int j4len, tlen;
3575     int nb, b;
3576
3577     /* All content of the new ci entry have already been filled correctly,
3578      * we only need to increase the count here (for non empty lists).
3579      */
3580     j4len = nbl->sci[nbl->nsci].cj4_ind_end - nbl->sci[nbl->nsci].cj4_ind_start;
3581     if (j4len > 0)
3582     {
3583         /* We can only have complete blocks of 4 j-entries in a list,
3584          * so round the count up before closing.
3585          */
3586         nbl->ncj4         = ((nbl->work->cj_ind + NBNXN_GPU_JGROUP_SIZE - 1) >> NBNXN_GPU_JGROUP_SIZE_2LOG);
3587         nbl->work->cj_ind = nbl->ncj4*NBNXN_GPU_JGROUP_SIZE;
3588
3589         nbl->nsci++;
3590
3591         if (nsp_max_av > 0)
3592         {
3593             /* Measure the size of the new entry and potentially split it */
3594             split_sci_entry(nbl, nsp_max_av, progBal, nc_bal, thread, nthread);
3595         }
3596     }
3597 }
3598
3599 /* Syncs the working array before adding another grid pair to the list */
3600 static void sync_work(nbnxn_pairlist_t *nbl)
3601 {
3602     if (!nbl->bSimple)
3603     {
3604         nbl->work->cj_ind   = nbl->ncj4*NBNXN_GPU_JGROUP_SIZE;
3605         nbl->work->cj4_init = nbl->ncj4;
3606     }
3607 }
3608
3609 /* Clears an nbnxn_pairlist_t data structure */
3610 static void clear_pairlist(nbnxn_pairlist_t *nbl)
3611 {
3612     nbl->nci           = 0;
3613     nbl->nsci          = 0;
3614     nbl->ncj           = 0;
3615     nbl->ncj4          = 0;
3616     nbl->nci_tot       = 0;
3617     nbl->nexcl         = 1;
3618
3619     nbl->work->ncj_noq = 0;
3620     nbl->work->ncj_hlj = 0;
3621 }
3622
3623 /* Sets a simple list i-cell bounding box, including PBC shift */
3624 static gmx_inline void set_icell_bb_simple(const nbnxn_bb_t *bb, int ci,
3625                                            real shx, real shy, real shz,
3626                                            nbnxn_bb_t *bb_ci)
3627 {
3628     bb_ci->lower[BB_X] = bb[ci].lower[BB_X] + shx;
3629     bb_ci->lower[BB_Y] = bb[ci].lower[BB_Y] + shy;
3630     bb_ci->lower[BB_Z] = bb[ci].lower[BB_Z] + shz;
3631     bb_ci->upper[BB_X] = bb[ci].upper[BB_X] + shx;
3632     bb_ci->upper[BB_Y] = bb[ci].upper[BB_Y] + shy;
3633     bb_ci->upper[BB_Z] = bb[ci].upper[BB_Z] + shz;
3634 }
3635
3636 #ifdef NBNXN_BBXXXX
3637 /* Sets a super-cell and sub cell bounding boxes, including PBC shift */
3638 static void set_icell_bbxxxx_supersub(const float *bb, int ci,
3639                                       real shx, real shy, real shz,
3640                                       float *bb_ci)
3641 {
3642     int ia, m, i;
3643
3644     ia = ci*(GPU_NSUBCELL>>STRIDE_PBB_2LOG)*NNBSBB_XXXX;
3645     for (m = 0; m < (GPU_NSUBCELL>>STRIDE_PBB_2LOG)*NNBSBB_XXXX; m += NNBSBB_XXXX)
3646     {
3647         for (i = 0; i < STRIDE_PBB; i++)
3648         {
3649             bb_ci[m+0*STRIDE_PBB+i] = bb[ia+m+0*STRIDE_PBB+i] + shx;
3650             bb_ci[m+1*STRIDE_PBB+i] = bb[ia+m+1*STRIDE_PBB+i] + shy;
3651             bb_ci[m+2*STRIDE_PBB+i] = bb[ia+m+2*STRIDE_PBB+i] + shz;
3652             bb_ci[m+3*STRIDE_PBB+i] = bb[ia+m+3*STRIDE_PBB+i] + shx;
3653             bb_ci[m+4*STRIDE_PBB+i] = bb[ia+m+4*STRIDE_PBB+i] + shy;
3654             bb_ci[m+5*STRIDE_PBB+i] = bb[ia+m+5*STRIDE_PBB+i] + shz;
3655         }
3656     }
3657 }
3658 #endif
3659
3660 /* Sets a super-cell and sub cell bounding boxes, including PBC shift */
3661 static void set_icell_bb_supersub(const nbnxn_bb_t *bb, int ci,
3662                                   real shx, real shy, real shz,
3663                                   nbnxn_bb_t *bb_ci)
3664 {
3665     int i;
3666
3667     for (i = 0; i < GPU_NSUBCELL; i++)
3668     {
3669         set_icell_bb_simple(bb, ci*GPU_NSUBCELL+i,
3670                             shx, shy, shz,
3671                             &bb_ci[i]);
3672     }
3673 }
3674
3675 /* Copies PBC shifted i-cell atom coordinates x,y,z to working array */
3676 static void icell_set_x_simple(int ci,
3677                                real shx, real shy, real shz,
3678                                int gmx_unused na_c,
3679                                int stride, const real *x,
3680                                nbnxn_list_work_t *work)
3681 {
3682     int  ia, i;
3683
3684     ia = ci*NBNXN_CPU_CLUSTER_I_SIZE;
3685
3686     for (i = 0; i < NBNXN_CPU_CLUSTER_I_SIZE; i++)
3687     {
3688         work->x_ci[i*STRIDE_XYZ+XX] = x[(ia+i)*stride+XX] + shx;
3689         work->x_ci[i*STRIDE_XYZ+YY] = x[(ia+i)*stride+YY] + shy;
3690         work->x_ci[i*STRIDE_XYZ+ZZ] = x[(ia+i)*stride+ZZ] + shz;
3691     }
3692 }
3693
3694 /* Copies PBC shifted super-cell atom coordinates x,y,z to working array */
3695 static void icell_set_x_supersub(int ci,
3696                                  real shx, real shy, real shz,
3697                                  int na_c,
3698                                  int stride, const real *x,
3699                                  nbnxn_list_work_t *work)
3700 {
3701     int  ia, i;
3702     real *x_ci;
3703
3704     x_ci = work->x_ci;
3705
3706     ia = ci*GPU_NSUBCELL*na_c;
3707     for (i = 0; i < GPU_NSUBCELL*na_c; i++)
3708     {
3709         x_ci[i*DIM + XX] = x[(ia+i)*stride + XX] + shx;
3710         x_ci[i*DIM + YY] = x[(ia+i)*stride + YY] + shy;
3711         x_ci[i*DIM + ZZ] = x[(ia+i)*stride + ZZ] + shz;
3712     }
3713 }
3714
3715 #ifdef NBNXN_SEARCH_BB_SIMD4
3716 /* Copies PBC shifted super-cell packed atom coordinates to working array */
3717 static void icell_set_x_supersub_simd4(int ci,
3718                                        real shx, real shy, real shz,
3719                                        int na_c,
3720                                        int stride, const real *x,
3721                                        nbnxn_list_work_t *work)
3722 {
3723     int  si, io, ia, i, j;
3724     real *x_ci;
3725
3726     x_ci = work->x_ci;
3727
3728     for (si = 0; si < GPU_NSUBCELL; si++)
3729     {
3730         for (i = 0; i < na_c; i += STRIDE_PBB)
3731         {
3732             io = si*na_c + i;
3733             ia = ci*GPU_NSUBCELL*na_c + io;
3734             for (j = 0; j < STRIDE_PBB; j++)
3735             {
3736                 x_ci[io*DIM + j + XX*STRIDE_PBB] = x[(ia+j)*stride+XX] + shx;
3737                 x_ci[io*DIM + j + YY*STRIDE_PBB] = x[(ia+j)*stride+YY] + shy;
3738                 x_ci[io*DIM + j + ZZ*STRIDE_PBB] = x[(ia+j)*stride+ZZ] + shz;
3739             }
3740         }
3741     }
3742 }
3743 #endif
3744
3745 static real nbnxn_rlist_inc_nonloc_fac = 0.6;
3746
3747 /* Due to the cluster size the effective pair-list is longer than
3748  * that of a simple atom pair-list. This function gives the extra distance.
3749  */
3750 real nbnxn_get_rlist_effective_inc(int cluster_size, real atom_density)
3751 {
3752     return ((0.5 + nbnxn_rlist_inc_nonloc_fac)*sqr(((cluster_size) - 1.0)/(cluster_size))*pow((cluster_size)/(atom_density), 1.0/3.0));
3753 }
3754
3755 /* Estimates the interaction volume^2 for non-local interactions */
3756 static real nonlocal_vol2(const gmx_domdec_zones_t *zones, rvec ls, real r)
3757 {
3758     int  z, d;
3759     real cl, ca, za;
3760     real vold_est;
3761     real vol2_est_tot;
3762
3763     vol2_est_tot = 0;
3764
3765     /* Here we simply add up the volumes of 1, 2 or 3 1D decomposition
3766      * not home interaction volume^2. As these volumes are not additive,
3767      * this is an overestimate, but it would only be significant in the limit
3768      * of small cells, where we anyhow need to split the lists into
3769      * as small parts as possible.
3770      */
3771
3772     for (z = 0; z < zones->n; z++)
3773     {
3774         if (zones->shift[z][XX] + zones->shift[z][YY] + zones->shift[z][ZZ] == 1)
3775         {
3776             cl = 0;
3777             ca = 1;
3778             za = 1;
3779             for (d = 0; d < DIM; d++)
3780             {
3781                 if (zones->shift[z][d] == 0)
3782                 {
3783                     cl += 0.5*ls[d];
3784                     ca *= ls[d];
3785                     za *= zones->size[z].x1[d] - zones->size[z].x0[d];
3786                 }
3787             }
3788
3789             /* 4 octants of a sphere */
3790             vold_est  = 0.25*M_PI*r*r*r*r;
3791             /* 4 quarter pie slices on the edges */
3792             vold_est += 4*cl*M_PI/6.0*r*r*r;
3793             /* One rectangular volume on a face */
3794             vold_est += ca*0.5*r*r;
3795
3796             vol2_est_tot += vold_est*za;
3797         }
3798     }
3799
3800     return vol2_est_tot;
3801 }
3802
3803 /* Estimates the average size of a full j-list for super/sub setup */
3804 static int get_nsubpair_max(const nbnxn_search_t nbs,
3805                             int                  iloc,
3806                             real                 rlist,
3807                             int                  min_ci_balanced)
3808 {
3809     const nbnxn_grid_t *grid;
3810     rvec ls;
3811     real xy_diag2, r_eff_sup, vol_est, nsp_est, nsp_est_nl;
3812     int  nsubpair_max;
3813
3814     grid = &nbs->grid[0];
3815
3816     ls[XX] = (grid->c1[XX] - grid->c0[XX])/(grid->ncx*GPU_NSUBCELL_X);
3817     ls[YY] = (grid->c1[YY] - grid->c0[YY])/(grid->ncy*GPU_NSUBCELL_Y);
3818     ls[ZZ] = (grid->c1[ZZ] - grid->c0[ZZ])*grid->ncx*grid->ncy/(grid->nc*GPU_NSUBCELL_Z);
3819
3820     /* The average squared length of the diagonal of a sub cell */
3821     xy_diag2 = ls[XX]*ls[XX] + ls[YY]*ls[YY] + ls[ZZ]*ls[ZZ];
3822
3823     /* The formulas below are a heuristic estimate of the average nsj per si*/
3824     r_eff_sup = rlist + nbnxn_rlist_inc_nonloc_fac*sqr((grid->na_c - 1.0)/grid->na_c)*sqrt(xy_diag2/3);
3825
3826     if (!nbs->DomDec || nbs->zones->n == 1)
3827     {
3828         nsp_est_nl = 0;
3829     }
3830     else
3831     {
3832         nsp_est_nl =
3833             sqr(grid->atom_density/grid->na_c)*
3834             nonlocal_vol2(nbs->zones, ls, r_eff_sup);
3835     }
3836
3837     if (LOCAL_I(iloc))
3838     {
3839         /* Sub-cell interacts with itself */
3840         vol_est  = ls[XX]*ls[YY]*ls[ZZ];
3841         /* 6/2 rectangular volume on the faces */
3842         vol_est += (ls[XX]*ls[YY] + ls[XX]*ls[ZZ] + ls[YY]*ls[ZZ])*r_eff_sup;
3843         /* 12/2 quarter pie slices on the edges */
3844         vol_est += 2*(ls[XX] + ls[YY] + ls[ZZ])*0.25*M_PI*sqr(r_eff_sup);
3845         /* 4 octants of a sphere */
3846         vol_est += 0.5*4.0/3.0*M_PI*pow(r_eff_sup, 3);
3847
3848         nsp_est = grid->nsubc_tot*vol_est*grid->atom_density/grid->na_c;
3849
3850         /* Subtract the non-local pair count */
3851         nsp_est -= nsp_est_nl;
3852
3853         if (debug)
3854         {
3855             fprintf(debug, "nsp_est local %5.1f non-local %5.1f\n",
3856                     nsp_est, nsp_est_nl);
3857         }
3858     }
3859     else
3860     {
3861         nsp_est = nsp_est_nl;
3862     }
3863
3864     if (min_ci_balanced <= 0 || grid->nc >= min_ci_balanced || grid->nc == 0)
3865     {
3866         /* We don't need to worry */
3867         nsubpair_max = -1;
3868     }
3869     else
3870     {
3871         /* Thus the (average) maximum j-list size should be as follows */
3872         nsubpair_max = max(1, (int)(nsp_est/min_ci_balanced+0.5));
3873
3874         /* Since the target value is a maximum (this avoids high outliers,
3875          * which lead to load imbalance), not average, we add half the
3876          * number of pairs in a cj4 block to get the average about right.
3877          */
3878         nsubpair_max += GPU_NSUBCELL*NBNXN_GPU_JGROUP_SIZE/2;
3879     }
3880
3881     if (debug)
3882     {
3883         fprintf(debug, "nbl nsp estimate %.1f, nsubpair_max %d\n",
3884                 nsp_est, nsubpair_max);
3885     }
3886
3887     return nsubpair_max;
3888 }
3889
3890 /* Debug list print function */
3891 static void print_nblist_ci_cj(FILE *fp, const nbnxn_pairlist_t *nbl)
3892 {
3893     int i, j;
3894
3895     for (i = 0; i < nbl->nci; i++)
3896     {
3897         fprintf(fp, "ci %4d  shift %2d  ncj %3d\n",
3898                 nbl->ci[i].ci, nbl->ci[i].shift,
3899                 nbl->ci[i].cj_ind_end - nbl->ci[i].cj_ind_start);
3900
3901         for (j = nbl->ci[i].cj_ind_start; j < nbl->ci[i].cj_ind_end; j++)
3902         {
3903             fprintf(fp, "  cj %5d  imask %x\n",
3904                     nbl->cj[j].cj,
3905                     nbl->cj[j].excl);
3906         }
3907     }
3908 }
3909
3910 /* Debug list print function */
3911 static void print_nblist_sci_cj(FILE *fp, const nbnxn_pairlist_t *nbl)
3912 {
3913     int i, j4, j, ncp, si;
3914
3915     for (i = 0; i < nbl->nsci; i++)
3916     {
3917         fprintf(fp, "ci %4d  shift %2d  ncj4 %2d\n",
3918                 nbl->sci[i].sci, nbl->sci[i].shift,
3919                 nbl->sci[i].cj4_ind_end - nbl->sci[i].cj4_ind_start);
3920
3921         ncp = 0;
3922         for (j4 = nbl->sci[i].cj4_ind_start; j4 < nbl->sci[i].cj4_ind_end; j4++)
3923         {
3924             for (j = 0; j < NBNXN_GPU_JGROUP_SIZE; j++)
3925             {
3926                 fprintf(fp, "  sj %5d  imask %x\n",
3927                         nbl->cj4[j4].cj[j],
3928                         nbl->cj4[j4].imei[0].imask);
3929                 for (si = 0; si < GPU_NSUBCELL; si++)
3930                 {
3931                     if (nbl->cj4[j4].imei[0].imask & (1U << (j*GPU_NSUBCELL + si)))
3932                     {
3933                         ncp++;
3934                     }
3935                 }
3936             }
3937         }
3938         fprintf(fp, "ci %4d  shift %2d  ncj4 %2d ncp %3d\n",
3939                 nbl->sci[i].sci, nbl->sci[i].shift,
3940                 nbl->sci[i].cj4_ind_end - nbl->sci[i].cj4_ind_start,
3941                 ncp);
3942     }
3943 }
3944
3945 /* Combine pair lists *nbl generated on multiple threads nblc */
3946 static void combine_nblists(int nnbl, nbnxn_pairlist_t **nbl,
3947                             nbnxn_pairlist_t *nblc)
3948 {
3949     int nsci, ncj4, nexcl;
3950     int n, i;
3951
3952     if (nblc->bSimple)
3953     {
3954         gmx_incons("combine_nblists does not support simple lists");
3955     }
3956
3957     nsci  = nblc->nsci;
3958     ncj4  = nblc->ncj4;
3959     nexcl = nblc->nexcl;
3960     for (i = 0; i < nnbl; i++)
3961     {
3962         nsci  += nbl[i]->nsci;
3963         ncj4  += nbl[i]->ncj4;
3964         nexcl += nbl[i]->nexcl;
3965     }
3966
3967     if (nsci > nblc->sci_nalloc)
3968     {
3969         nb_realloc_sci(nblc, nsci);
3970     }
3971     if (ncj4 > nblc->cj4_nalloc)
3972     {
3973         nblc->cj4_nalloc = over_alloc_small(ncj4);
3974         nbnxn_realloc_void((void **)&nblc->cj4,
3975                            nblc->ncj4*sizeof(*nblc->cj4),
3976                            nblc->cj4_nalloc*sizeof(*nblc->cj4),
3977                            nblc->alloc, nblc->free);
3978     }
3979     if (nexcl > nblc->excl_nalloc)
3980     {
3981         nblc->excl_nalloc = over_alloc_small(nexcl);
3982         nbnxn_realloc_void((void **)&nblc->excl,
3983                            nblc->nexcl*sizeof(*nblc->excl),
3984                            nblc->excl_nalloc*sizeof(*nblc->excl),
3985                            nblc->alloc, nblc->free);
3986     }
3987
3988     /* Each thread should copy its own data to the combined arrays,
3989      * as otherwise data will go back and forth between different caches.
3990      */
3991 #pragma omp parallel for num_threads(gmx_omp_nthreads_get(emntPairsearch)) schedule(static)
3992     for (n = 0; n < nnbl; n++)
3993     {
3994         int sci_offset;
3995         int cj4_offset;
3996         int ci_offset;
3997         int excl_offset;
3998         int i, j4;
3999         const nbnxn_pairlist_t *nbli;
4000
4001         /* Determine the offset in the combined data for our thread */
4002         sci_offset  = nblc->nsci;
4003         cj4_offset  = nblc->ncj4;
4004         ci_offset   = nblc->nci_tot;
4005         excl_offset = nblc->nexcl;
4006
4007         for (i = 0; i < n; i++)
4008         {
4009             sci_offset  += nbl[i]->nsci;
4010             cj4_offset  += nbl[i]->ncj4;
4011             ci_offset   += nbl[i]->nci_tot;
4012             excl_offset += nbl[i]->nexcl;
4013         }
4014
4015         nbli = nbl[n];
4016
4017         for (i = 0; i < nbli->nsci; i++)
4018         {
4019             nblc->sci[sci_offset+i]                = nbli->sci[i];
4020             nblc->sci[sci_offset+i].cj4_ind_start += cj4_offset;
4021             nblc->sci[sci_offset+i].cj4_ind_end   += cj4_offset;
4022         }
4023
4024         for (j4 = 0; j4 < nbli->ncj4; j4++)
4025         {
4026             nblc->cj4[cj4_offset+j4]                   = nbli->cj4[j4];
4027             nblc->cj4[cj4_offset+j4].imei[0].excl_ind += excl_offset;
4028             nblc->cj4[cj4_offset+j4].imei[1].excl_ind += excl_offset;
4029         }
4030
4031         for (j4 = 0; j4 < nbli->nexcl; j4++)
4032         {
4033             nblc->excl[excl_offset+j4] = nbli->excl[j4];
4034         }
4035     }
4036
4037     for (n = 0; n < nnbl; n++)
4038     {
4039         nblc->nsci    += nbl[n]->nsci;
4040         nblc->ncj4    += nbl[n]->ncj4;
4041         nblc->nci_tot += nbl[n]->nci_tot;
4042         nblc->nexcl   += nbl[n]->nexcl;
4043     }
4044 }
4045
4046 /* Returns the next ci to be processes by our thread */
4047 static gmx_bool next_ci(const nbnxn_grid_t *grid,
4048                         int conv,
4049                         int nth, int ci_block,
4050                         int *ci_x, int *ci_y,
4051                         int *ci_b, int *ci)
4052 {
4053     (*ci_b)++;
4054     (*ci)++;
4055
4056     if (*ci_b == ci_block)
4057     {
4058         /* Jump to the next block assigned to this task */
4059         *ci   += (nth - 1)*ci_block;
4060         *ci_b  = 0;
4061     }
4062
4063     if (*ci >= grid->nc*conv)
4064     {
4065         return FALSE;
4066     }
4067
4068     while (*ci >= grid->cxy_ind[*ci_x*grid->ncy + *ci_y + 1]*conv)
4069     {
4070         *ci_y += 1;
4071         if (*ci_y == grid->ncy)
4072         {
4073             *ci_x += 1;
4074             *ci_y  = 0;
4075         }
4076     }
4077
4078     return TRUE;
4079 }
4080
4081 /* Returns the distance^2 for which we put cell pairs in the list
4082  * without checking atom pair distances. This is usually < rlist^2.
4083  */
4084 static float boundingbox_only_distance2(const nbnxn_grid_t *gridi,
4085                                         const nbnxn_grid_t *gridj,
4086                                         real                rlist,
4087                                         gmx_bool            simple)
4088 {
4089     /* If the distance between two sub-cell bounding boxes is less
4090      * than this distance, do not check the distance between
4091      * all particle pairs in the sub-cell, since then it is likely
4092      * that the box pair has atom pairs within the cut-off.
4093      * We use the nblist cut-off minus 0.5 times the average x/y diagonal
4094      * spacing of the sub-cells. Around 40% of the checked pairs are pruned.
4095      * Using more than 0.5 gains at most 0.5%.
4096      * If forces are calculated more than twice, the performance gain
4097      * in the force calculation outweighs the cost of checking.
4098      * Note that with subcell lists, the atom-pair distance check
4099      * is only performed when only 1 out of 8 sub-cells in within range,
4100      * this is because the GPU is much faster than the cpu.
4101      */
4102     real bbx, bby;
4103     real rbb2;
4104
4105     bbx = 0.5*(gridi->sx + gridj->sx);
4106     bby = 0.5*(gridi->sy + gridj->sy);
4107     if (!simple)
4108     {
4109         bbx /= GPU_NSUBCELL_X;
4110         bby /= GPU_NSUBCELL_Y;
4111     }
4112
4113     rbb2 = sqr(max(0, rlist - 0.5*sqrt(bbx*bbx + bby*bby)));
4114
4115 #ifndef GMX_DOUBLE
4116     return rbb2;
4117 #else
4118     return (float)((1+GMX_FLOAT_EPS)*rbb2);
4119 #endif
4120 }
4121
4122 static int get_ci_block_size(const nbnxn_grid_t *gridi,
4123                              gmx_bool bDomDec, int nth)
4124 {
4125     const int ci_block_enum      = 5;
4126     const int ci_block_denom     = 11;
4127     const int ci_block_min_atoms = 16;
4128     int ci_block;
4129
4130     /* Here we decide how to distribute the blocks over the threads.
4131      * We use prime numbers to try to avoid that the grid size becomes
4132      * a multiple of the number of threads, which would lead to some
4133      * threads getting "inner" pairs and others getting boundary pairs,
4134      * which in turns will lead to load imbalance between threads.
4135      * Set the block size as 5/11/ntask times the average number of cells
4136      * in a y,z slab. This should ensure a quite uniform distribution
4137      * of the grid parts of the different thread along all three grid
4138      * zone boundaries with 3D domain decomposition. At the same time
4139      * the blocks will not become too small.
4140      */
4141     ci_block = (gridi->nc*ci_block_enum)/(ci_block_denom*gridi->ncx*nth);
4142
4143     /* Ensure the blocks are not too small: avoids cache invalidation */
4144     if (ci_block*gridi->na_sc < ci_block_min_atoms)
4145     {
4146         ci_block = (ci_block_min_atoms + gridi->na_sc - 1)/gridi->na_sc;
4147     }
4148
4149     /* Without domain decomposition
4150      * or with less than 3 blocks per task, divide in nth blocks.
4151      */
4152     if (!bDomDec || ci_block*3*nth > gridi->nc)
4153     {
4154         ci_block = (gridi->nc + nth - 1)/nth;
4155     }
4156
4157     return ci_block;
4158 }
4159
4160 /* Generates the part of pair-list nbl assigned to our thread */
4161 static void nbnxn_make_pairlist_part(const nbnxn_search_t nbs,
4162                                      const nbnxn_grid_t *gridi,
4163                                      const nbnxn_grid_t *gridj,
4164                                      nbnxn_search_work_t *work,
4165                                      const nbnxn_atomdata_t *nbat,
4166                                      const t_blocka *excl,
4167                                      real rlist,
4168                                      int nb_kernel_type,
4169                                      int ci_block,
4170                                      gmx_bool bFBufferFlag,
4171                                      int nsubpair_max,
4172                                      gmx_bool progBal,
4173                                      int min_ci_balanced,
4174                                      int th, int nth,
4175                                      nbnxn_pairlist_t *nbl)
4176 {
4177     int  na_cj_2log;
4178     matrix box;
4179     real rl2;
4180     float rbb2;
4181     int  d;
4182     int  ci_b, ci, ci_x, ci_y, ci_xy, cj;
4183     ivec shp;
4184     int  tx, ty, tz;
4185     int  shift;
4186     gmx_bool bMakeList;
4187     real shx, shy, shz;
4188     int  conv_i, cell0_i;
4189     const nbnxn_bb_t *bb_i = NULL;
4190 #ifdef NBNXN_BBXXXX
4191     const float *pbb_i = NULL;
4192 #endif
4193     const float *bbcz_i, *bbcz_j;
4194     const int *flags_i;
4195     real bx0, bx1, by0, by1, bz0, bz1;
4196     real bz1_frac;
4197     real d2cx, d2z, d2z_cx, d2z_cy, d2zx, d2zxy, d2xy;
4198     int  cxf, cxl, cyf, cyf_x, cyl;
4199     int  cx, cy;
4200     int  c0, c1, cs, cf, cl;
4201     int  ndistc;
4202     int  ncpcheck;
4203     int  gridi_flag_shift = 0, gridj_flag_shift = 0;
4204     unsigned *gridj_flag  = NULL;
4205     int  ncj_old_i, ncj_old_j;
4206
4207     nbs_cycle_start(&work->cc[enbsCCsearch]);
4208
4209     if (gridj->bSimple != nbl->bSimple)
4210     {
4211         gmx_incons("Grid incompatible with pair-list");
4212     }
4213
4214     sync_work(nbl);
4215     nbl->na_sc = gridj->na_sc;
4216     nbl->na_ci = gridj->na_c;
4217     nbl->na_cj = nbnxn_kernel_to_cj_size(nb_kernel_type);
4218     na_cj_2log = get_2log(nbl->na_cj);
4219
4220     nbl->rlist  = rlist;
4221
4222     if (bFBufferFlag)
4223     {
4224         /* Determine conversion of clusters to flag blocks */
4225         gridi_flag_shift = 0;
4226         while ((nbl->na_ci<<gridi_flag_shift) < NBNXN_BUFFERFLAG_SIZE)
4227         {
4228             gridi_flag_shift++;
4229         }
4230         gridj_flag_shift = 0;
4231         while ((nbl->na_cj<<gridj_flag_shift) < NBNXN_BUFFERFLAG_SIZE)
4232         {
4233             gridj_flag_shift++;
4234         }
4235
4236         gridj_flag = work->buffer_flags.flag;
4237     }
4238
4239     copy_mat(nbs->box, box);
4240
4241     rl2 = nbl->rlist*nbl->rlist;
4242
4243     rbb2 = boundingbox_only_distance2(gridi, gridj, nbl->rlist, nbl->bSimple);
4244
4245     if (debug)
4246     {
4247         fprintf(debug, "nbl bounding box only distance %f\n", sqrt(rbb2));
4248     }
4249
4250     /* Set the shift range */
4251     for (d = 0; d < DIM; d++)
4252     {
4253         /* Check if we need periodicity shifts.
4254          * Without PBC or with domain decomposition we don't need them.
4255          */
4256         if (d >= ePBC2npbcdim(nbs->ePBC) || nbs->dd_dim[d])
4257         {
4258             shp[d] = 0;
4259         }
4260         else
4261         {
4262             if (d == XX &&
4263                 box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < sqrt(rl2))
4264             {
4265                 shp[d] = 2;
4266             }
4267             else
4268             {
4269                 shp[d] = 1;
4270             }
4271         }
4272     }
4273
4274     if (nbl->bSimple && !gridi->bSimple)
4275     {
4276         conv_i  = gridi->na_sc/gridj->na_sc;
4277         bb_i    = gridi->bb_simple;
4278         bbcz_i  = gridi->bbcz_simple;
4279         flags_i = gridi->flags_simple;
4280     }
4281     else
4282     {
4283         conv_i  = 1;
4284 #ifdef NBNXN_BBXXXX
4285         if (gridi->bSimple)
4286         {
4287             bb_i  = gridi->bb;
4288         }
4289         else
4290         {
4291             pbb_i = gridi->pbb;
4292         }
4293 #else
4294         /* We use the normal bounding box format for both grid types */
4295         bb_i  = gridi->bb;
4296 #endif
4297         bbcz_i  = gridi->bbcz;
4298         flags_i = gridi->flags;
4299     }
4300     cell0_i = gridi->cell0*conv_i;
4301
4302     bbcz_j = gridj->bbcz;
4303
4304     if (conv_i != 1)
4305     {
4306         /* Blocks of the conversion factor - 1 give a large repeat count
4307          * combined with a small block size. This should result in good
4308          * load balancing for both small and large domains.
4309          */
4310         ci_block = conv_i - 1;
4311     }
4312     if (debug)
4313     {
4314         fprintf(debug, "nbl nc_i %d col.av. %.1f ci_block %d\n",
4315                 gridi->nc, gridi->nc/(double)(gridi->ncx*gridi->ncy), ci_block);
4316     }
4317
4318     ndistc   = 0;
4319     ncpcheck = 0;
4320
4321     /* Initially ci_b and ci to 1 before where we want them to start,
4322      * as they will both be incremented in next_ci.
4323      */
4324     ci_b = -1;
4325     ci   = th*ci_block - 1;
4326     ci_x = 0;
4327     ci_y = 0;
4328     while (next_ci(gridi, conv_i, nth, ci_block, &ci_x, &ci_y, &ci_b, &ci))
4329     {
4330         if (nbl->bSimple && flags_i[ci] == 0)
4331         {
4332             continue;
4333         }
4334
4335         ncj_old_i = nbl->ncj;
4336
4337         d2cx = 0;
4338         if (gridj != gridi && shp[XX] == 0)
4339         {
4340             if (nbl->bSimple)
4341             {
4342                 bx1 = bb_i[ci].upper[BB_X];
4343             }
4344             else
4345             {
4346                 bx1 = gridi->c0[XX] + (ci_x+1)*gridi->sx;
4347             }
4348             if (bx1 < gridj->c0[XX])
4349             {
4350                 d2cx = sqr(gridj->c0[XX] - bx1);
4351
4352                 if (d2cx >= rl2)
4353                 {
4354                     continue;
4355                 }
4356             }
4357         }
4358
4359         ci_xy = ci_x*gridi->ncy + ci_y;
4360
4361         /* Loop over shift vectors in three dimensions */
4362         for (tz = -shp[ZZ]; tz <= shp[ZZ]; tz++)
4363         {
4364             shz = tz*box[ZZ][ZZ];
4365
4366             bz0 = bbcz_i[ci*NNBSBB_D  ] + shz;
4367             bz1 = bbcz_i[ci*NNBSBB_D+1] + shz;
4368
4369             if (tz == 0)
4370             {
4371                 d2z = 0;
4372             }
4373             else if (tz < 0)
4374             {
4375                 d2z = sqr(bz1);
4376             }
4377             else
4378             {
4379                 d2z = sqr(bz0 - box[ZZ][ZZ]);
4380             }
4381
4382             d2z_cx = d2z + d2cx;
4383
4384             if (d2z_cx >= rl2)
4385             {
4386                 continue;
4387             }
4388
4389             bz1_frac =
4390                 bz1/((real)(gridi->cxy_ind[ci_xy+1] - gridi->cxy_ind[ci_xy]));
4391             if (bz1_frac < 0)
4392             {
4393                 bz1_frac = 0;
4394             }
4395             /* The check with bz1_frac close to or larger than 1 comes later */
4396
4397             for (ty = -shp[YY]; ty <= shp[YY]; ty++)
4398             {
4399                 shy = ty*box[YY][YY] + tz*box[ZZ][YY];
4400
4401                 if (nbl->bSimple)
4402                 {
4403                     by0 = bb_i[ci].lower[BB_Y] + shy;
4404                     by1 = bb_i[ci].upper[BB_Y] + shy;
4405                 }
4406                 else
4407                 {
4408                     by0 = gridi->c0[YY] + (ci_y  )*gridi->sy + shy;
4409                     by1 = gridi->c0[YY] + (ci_y+1)*gridi->sy + shy;
4410                 }
4411
4412                 get_cell_range(by0, by1,
4413                                gridj->ncy, gridj->c0[YY], gridj->sy, gridj->inv_sy,
4414                                d2z_cx, rl2,
4415                                &cyf, &cyl);
4416
4417                 if (cyf > cyl)
4418                 {
4419                     continue;
4420                 }
4421
4422                 d2z_cy = d2z;
4423                 if (by1 < gridj->c0[YY])
4424                 {
4425                     d2z_cy += sqr(gridj->c0[YY] - by1);
4426                 }
4427                 else if (by0 > gridj->c1[YY])
4428                 {
4429                     d2z_cy += sqr(by0 - gridj->c1[YY]);
4430                 }
4431
4432                 for (tx = -shp[XX]; tx <= shp[XX]; tx++)
4433                 {
4434                     shift = XYZ2IS(tx, ty, tz);
4435
4436 #ifdef NBNXN_SHIFT_BACKWARD
4437                     if (gridi == gridj && shift > CENTRAL)
4438                     {
4439                         continue;
4440                     }
4441 #endif
4442
4443                     shx = tx*box[XX][XX] + ty*box[YY][XX] + tz*box[ZZ][XX];
4444
4445                     if (nbl->bSimple)
4446                     {
4447                         bx0 = bb_i[ci].lower[BB_X] + shx;
4448                         bx1 = bb_i[ci].upper[BB_X] + shx;
4449                     }
4450                     else
4451                     {
4452                         bx0 = gridi->c0[XX] + (ci_x  )*gridi->sx + shx;
4453                         bx1 = gridi->c0[XX] + (ci_x+1)*gridi->sx + shx;
4454                     }
4455
4456                     get_cell_range(bx0, bx1,
4457                                    gridj->ncx, gridj->c0[XX], gridj->sx, gridj->inv_sx,
4458                                    d2z_cy, rl2,
4459                                    &cxf, &cxl);
4460
4461                     if (cxf > cxl)
4462                     {
4463                         continue;
4464                     }
4465
4466                     if (nbl->bSimple)
4467                     {
4468                         new_ci_entry(nbl, cell0_i+ci, shift, flags_i[ci]);
4469                     }
4470                     else
4471                     {
4472                         new_sci_entry(nbl, cell0_i+ci, shift);
4473                     }
4474
4475 #ifndef NBNXN_SHIFT_BACKWARD
4476                     if (cxf < ci_x)
4477 #else
4478                     if (shift == CENTRAL && gridi == gridj &&
4479                         cxf < ci_x)
4480 #endif
4481                     {
4482                         /* Leave the pairs with i > j.
4483                          * x is the major index, so skip half of it.
4484                          */
4485                         cxf = ci_x;
4486                     }
4487
4488                     if (nbl->bSimple)
4489                     {
4490                         set_icell_bb_simple(bb_i, ci, shx, shy, shz,
4491                                             nbl->work->bb_ci);
4492                     }
4493                     else
4494                     {
4495 #ifdef NBNXN_BBXXXX
4496                         set_icell_bbxxxx_supersub(pbb_i, ci, shx, shy, shz,
4497                                                   nbl->work->pbb_ci);
4498 #else
4499                         set_icell_bb_supersub(bb_i, ci, shx, shy, shz,
4500                                               nbl->work->bb_ci);
4501 #endif
4502                     }
4503
4504                     nbs->icell_set_x(cell0_i+ci, shx, shy, shz,
4505                                      gridi->na_c, nbat->xstride, nbat->x,
4506                                      nbl->work);
4507
4508                     for (cx = cxf; cx <= cxl; cx++)
4509                     {
4510                         d2zx = d2z;
4511                         if (gridj->c0[XX] + cx*gridj->sx > bx1)
4512                         {
4513                             d2zx += sqr(gridj->c0[XX] + cx*gridj->sx - bx1);
4514                         }
4515                         else if (gridj->c0[XX] + (cx+1)*gridj->sx < bx0)
4516                         {
4517                             d2zx += sqr(gridj->c0[XX] + (cx+1)*gridj->sx - bx0);
4518                         }
4519
4520 #ifndef NBNXN_SHIFT_BACKWARD
4521                         if (gridi == gridj &&
4522                             cx == 0 && cyf < ci_y)
4523 #else
4524                         if (gridi == gridj &&
4525                             cx == 0 && shift == CENTRAL && cyf < ci_y)
4526 #endif
4527                         {
4528                             /* Leave the pairs with i > j.
4529                              * Skip half of y when i and j have the same x.
4530                              */
4531                             cyf_x = ci_y;
4532                         }
4533                         else
4534                         {
4535                             cyf_x = cyf;
4536                         }
4537
4538                         for (cy = cyf_x; cy <= cyl; cy++)
4539                         {
4540                             c0 = gridj->cxy_ind[cx*gridj->ncy+cy];
4541                             c1 = gridj->cxy_ind[cx*gridj->ncy+cy+1];
4542 #ifdef NBNXN_SHIFT_BACKWARD
4543                             if (gridi == gridj &&
4544                                 shift == CENTRAL && c0 < ci)
4545                             {
4546                                 c0 = ci;
4547                             }
4548 #endif
4549
4550                             d2zxy = d2zx;
4551                             if (gridj->c0[YY] + cy*gridj->sy > by1)
4552                             {
4553                                 d2zxy += sqr(gridj->c0[YY] + cy*gridj->sy - by1);
4554                             }
4555                             else if (gridj->c0[YY] + (cy+1)*gridj->sy < by0)
4556                             {
4557                                 d2zxy += sqr(gridj->c0[YY] + (cy+1)*gridj->sy - by0);
4558                             }
4559                             if (c1 > c0 && d2zxy < rl2)
4560                             {
4561                                 cs = c0 + (int)(bz1_frac*(c1 - c0));
4562                                 if (cs >= c1)
4563                                 {
4564                                     cs = c1 - 1;
4565                                 }
4566
4567                                 d2xy = d2zxy - d2z;
4568
4569                                 /* Find the lowest cell that can possibly
4570                                  * be within range.
4571                                  */
4572                                 cf = cs;
4573                                 while (cf > c0 &&
4574                                        (bbcz_j[cf*NNBSBB_D+1] >= bz0 ||
4575                                         d2xy + sqr(bbcz_j[cf*NNBSBB_D+1] - bz0) < rl2))
4576                                 {
4577                                     cf--;
4578                                 }
4579
4580                                 /* Find the highest cell that can possibly
4581                                  * be within range.
4582                                  */
4583                                 cl = cs;
4584                                 while (cl < c1-1 &&
4585                                        (bbcz_j[cl*NNBSBB_D] <= bz1 ||
4586                                         d2xy + sqr(bbcz_j[cl*NNBSBB_D] - bz1) < rl2))
4587                                 {
4588                                     cl++;
4589                                 }
4590
4591 #ifdef NBNXN_REFCODE
4592                                 {
4593                                     /* Simple reference code, for debugging,
4594                                      * overrides the more complex code above.
4595                                      */
4596                                     int k;
4597                                     cf = c1;
4598                                     cl = -1;
4599                                     for (k = c0; k < c1; k++)
4600                                     {
4601                                         if (box_dist2(bx0, bx1, by0, by1, bz0, bz1, bb+k) < rl2 &&
4602                                             k < cf)
4603                                         {
4604                                             cf = k;
4605                                         }
4606                                         if (box_dist2(bx0, bx1, by0, by1, bz0, bz1, bb+k) < rl2 &&
4607                                             k > cl)
4608                                         {
4609                                             cl = k;
4610                                         }
4611                                     }
4612                                 }
4613 #endif
4614
4615                                 if (gridi == gridj)
4616                                 {
4617                                     /* We want each atom/cell pair only once,
4618                                      * only use cj >= ci.
4619                                      */
4620 #ifndef NBNXN_SHIFT_BACKWARD
4621                                     cf = max(cf, ci);
4622 #else
4623                                     if (shift == CENTRAL)
4624                                     {
4625                                         cf = max(cf, ci);
4626                                     }
4627 #endif
4628                                 }
4629
4630                                 if (cf <= cl)
4631                                 {
4632                                     /* For f buffer flags with simple lists */
4633                                     ncj_old_j = nbl->ncj;
4634
4635                                     switch (nb_kernel_type)
4636                                     {
4637                                         case nbnxnk4x4_PlainC:
4638                                             check_subcell_list_space_simple(nbl, cl-cf+1);
4639
4640                                             make_cluster_list_simple(gridj,
4641                                                                      nbl, ci, cf, cl,
4642                                                                      (gridi == gridj && shift == CENTRAL),
4643                                                                      nbat->x,
4644                                                                      rl2, rbb2,
4645                                                                      &ndistc);
4646                                             break;
4647 #ifdef GMX_NBNXN_SIMD_4XN
4648                                         case nbnxnk4xN_SIMD_4xN:
4649                                             check_subcell_list_space_simple(nbl, ci_to_cj(na_cj_2log, cl-cf)+2);
4650                                             make_cluster_list_simd_4xn(gridj,
4651                                                                        nbl, ci, cf, cl,
4652                                                                        (gridi == gridj && shift == CENTRAL),
4653                                                                        nbat->x,
4654                                                                        rl2, rbb2,
4655                                                                        &ndistc);
4656                                             break;
4657 #endif
4658 #ifdef GMX_NBNXN_SIMD_2XNN
4659                                         case nbnxnk4xN_SIMD_2xNN:
4660                                             check_subcell_list_space_simple(nbl, ci_to_cj(na_cj_2log, cl-cf)+2);
4661                                             make_cluster_list_simd_2xnn(gridj,
4662                                                                         nbl, ci, cf, cl,
4663                                                                         (gridi == gridj && shift == CENTRAL),
4664                                                                         nbat->x,
4665                                                                         rl2, rbb2,
4666                                                                         &ndistc);
4667                                             break;
4668 #endif
4669                                         case nbnxnk8x8x8_PlainC:
4670                                         case nbnxnk8x8x8_CUDA:
4671                                             check_subcell_list_space_supersub(nbl, cl-cf+1);
4672                                             for (cj = cf; cj <= cl; cj++)
4673                                             {
4674                                                 make_cluster_list_supersub(gridi, gridj,
4675                                                                            nbl, ci, cj,
4676                                                                            (gridi == gridj && shift == CENTRAL && ci == cj),
4677                                                                            nbat->xstride, nbat->x,
4678                                                                            rl2, rbb2,
4679                                                                            &ndistc);
4680                                             }
4681                                             break;
4682                                     }
4683                                     ncpcheck += cl - cf + 1;
4684
4685                                     if (bFBufferFlag && nbl->ncj > ncj_old_j)
4686                                     {
4687                                         int cbf, cbl, cb;
4688
4689                                         cbf = nbl->cj[ncj_old_j].cj >> gridj_flag_shift;
4690                                         cbl = nbl->cj[nbl->ncj-1].cj >> gridj_flag_shift;
4691                                         for (cb = cbf; cb <= cbl; cb++)
4692                                         {
4693                                             gridj_flag[cb] = 1U<<th;
4694                                         }
4695                                     }
4696                                 }
4697                             }
4698                         }
4699                     }
4700
4701                     /* Set the exclusions for this ci list */
4702                     if (nbl->bSimple)
4703                     {
4704                         set_ci_top_excls(nbs,
4705                                          nbl,
4706                                          shift == CENTRAL && gridi == gridj,
4707                                          gridj->na_c_2log,
4708                                          na_cj_2log,
4709                                          &(nbl->ci[nbl->nci]),
4710                                          excl);
4711                     }
4712                     else
4713                     {
4714                         set_sci_top_excls(nbs,
4715                                           nbl,
4716                                           shift == CENTRAL && gridi == gridj,
4717                                           gridj->na_c_2log,
4718                                           &(nbl->sci[nbl->nsci]),
4719                                           excl);
4720                     }
4721
4722                     /* Close this ci list */
4723                     if (nbl->bSimple)
4724                     {
4725                         close_ci_entry_simple(nbl);
4726                     }
4727                     else
4728                     {
4729                         close_ci_entry_supersub(nbl,
4730                                                 nsubpair_max,
4731                                                 progBal, min_ci_balanced,
4732                                                 th, nth);
4733                     }
4734                 }
4735             }
4736         }
4737
4738         if (bFBufferFlag && nbl->ncj > ncj_old_i)
4739         {
4740             work->buffer_flags.flag[(gridi->cell0+ci)>>gridi_flag_shift] = 1U<<th;
4741         }
4742     }
4743
4744     work->ndistc = ndistc;
4745
4746     nbs_cycle_stop(&work->cc[enbsCCsearch]);
4747
4748     if (debug)
4749     {
4750         fprintf(debug, "number of distance checks %d\n", ndistc);
4751         fprintf(debug, "ncpcheck %s %d\n", gridi == gridj ? "local" : "non-local",
4752                 ncpcheck);
4753
4754         if (nbl->bSimple)
4755         {
4756             print_nblist_statistics_simple(debug, nbl, nbs, rlist);
4757         }
4758         else
4759         {
4760             print_nblist_statistics_supersub(debug, nbl, nbs, rlist);
4761         }
4762
4763     }
4764 }
4765
4766 static void reduce_buffer_flags(const nbnxn_search_t        nbs,
4767                                 int                         nsrc,
4768                                 const nbnxn_buffer_flags_t *dest)
4769 {
4770     int s, b;
4771     const unsigned *flag;
4772
4773     for (s = 0; s < nsrc; s++)
4774     {
4775         flag = nbs->work[s].buffer_flags.flag;
4776
4777         for (b = 0; b < dest->nflag; b++)
4778         {
4779             dest->flag[b] |= flag[b];
4780         }
4781     }
4782 }
4783
4784 static void print_reduction_cost(const nbnxn_buffer_flags_t *flags, int nout)
4785 {
4786     int nelem, nkeep, ncopy, nred, b, c, out;
4787
4788     nelem = 0;
4789     nkeep = 0;
4790     ncopy = 0;
4791     nred  = 0;
4792     for (b = 0; b < flags->nflag; b++)
4793     {
4794         if (flags->flag[b] == 1)
4795         {
4796             /* Only flag 0 is set, no copy of reduction required */
4797             nelem++;
4798             nkeep++;
4799         }
4800         else if (flags->flag[b] > 0)
4801         {
4802             c = 0;
4803             for (out = 0; out < nout; out++)
4804             {
4805                 if (flags->flag[b] & (1U<<out))
4806                 {
4807                     c++;
4808                 }
4809             }
4810             nelem += c;
4811             if (c == 1)
4812             {
4813                 ncopy++;
4814             }
4815             else
4816             {
4817                 nred += c;
4818             }
4819         }
4820     }
4821
4822     fprintf(debug, "nbnxn reduction: #flag %d #list %d elem %4.2f, keep %4.2f copy %4.2f red %4.2f\n",
4823             flags->nflag, nout,
4824             nelem/(double)(flags->nflag),
4825             nkeep/(double)(flags->nflag),
4826             ncopy/(double)(flags->nflag),
4827             nred/(double)(flags->nflag));
4828 }
4829
4830 /* Perform a count (linear) sort to sort the smaller lists to the end.
4831  * This avoids load imbalance on the GPU, as large lists will be
4832  * scheduled and executed first and the smaller lists later.
4833  * Load balancing between multi-processors only happens at the end
4834  * and there smaller lists lead to more effective load balancing.
4835  * The sorting is done on the cj4 count, not on the actual pair counts.
4836  * Not only does this make the sort faster, but it also results in
4837  * better load balancing than using a list sorted on exact load.
4838  * This function swaps the pointer in the pair list to avoid a copy operation.
4839  */
4840 static void sort_sci(nbnxn_pairlist_t *nbl)
4841 {
4842     nbnxn_list_work_t *work;
4843     int                m, i, s, s0, s1;
4844     nbnxn_sci_t       *sci_sort;
4845
4846     if (nbl->ncj4 <= nbl->nsci)
4847     {
4848         /* nsci = 0 or all sci have size 1, sorting won't change the order */
4849         return;
4850     }
4851
4852     work = nbl->work;
4853
4854     /* We will distinguish differences up to double the average */
4855     m = (2*nbl->ncj4)/nbl->nsci;
4856
4857     if (m + 1 > work->sort_nalloc)
4858     {
4859         work->sort_nalloc = over_alloc_large(m + 1);
4860         srenew(work->sort, work->sort_nalloc);
4861     }
4862
4863     if (work->sci_sort_nalloc != nbl->sci_nalloc)
4864     {
4865         work->sci_sort_nalloc = nbl->sci_nalloc;
4866         nbnxn_realloc_void((void **)&work->sci_sort,
4867                            0,
4868                            work->sci_sort_nalloc*sizeof(*work->sci_sort),
4869                            nbl->alloc, nbl->free);
4870     }
4871
4872     /* Count the entries of each size */
4873     for (i = 0; i <= m; i++)
4874     {
4875         work->sort[i] = 0;
4876     }
4877     for (s = 0; s < nbl->nsci; s++)
4878     {
4879         i = min(m, nbl->sci[s].cj4_ind_end - nbl->sci[s].cj4_ind_start);
4880         work->sort[i]++;
4881     }
4882     /* Calculate the offset for each count */
4883     s0            = work->sort[m];
4884     work->sort[m] = 0;
4885     for (i = m - 1; i >= 0; i--)
4886     {
4887         s1            = work->sort[i];
4888         work->sort[i] = work->sort[i + 1] + s0;
4889         s0            = s1;
4890     }
4891
4892     /* Sort entries directly into place */
4893     sci_sort = work->sci_sort;
4894     for (s = 0; s < nbl->nsci; s++)
4895     {
4896         i = min(m, nbl->sci[s].cj4_ind_end - nbl->sci[s].cj4_ind_start);
4897         sci_sort[work->sort[i]++] = nbl->sci[s];
4898     }
4899
4900     /* Swap the sci pointers so we use the new, sorted list */
4901     work->sci_sort = nbl->sci;
4902     nbl->sci       = sci_sort;
4903 }
4904
4905 /* Make a local or non-local pair-list, depending on iloc */
4906 void nbnxn_make_pairlist(const nbnxn_search_t  nbs,
4907                          nbnxn_atomdata_t     *nbat,
4908                          const t_blocka       *excl,
4909                          real                  rlist,
4910                          int                   min_ci_balanced,
4911                          nbnxn_pairlist_set_t *nbl_list,
4912                          int                   iloc,
4913                          int                   nb_kernel_type,
4914                          t_nrnb               *nrnb)
4915 {
4916     nbnxn_grid_t *gridi, *gridj;
4917     gmx_bool bGPUCPU;
4918     int nzi, zi, zj0, zj1, zj;
4919     int nsubpair_max;
4920     int th;
4921     int nnbl;
4922     nbnxn_pairlist_t **nbl;
4923     int ci_block;
4924     gmx_bool CombineNBLists;
4925     gmx_bool progBal;
4926     int np_tot, np_noq, np_hlj, nap;
4927
4928     /* Check if we are running hybrid GPU + CPU nbnxn mode */
4929     bGPUCPU = (!nbs->grid[0].bSimple && nbl_list->bSimple);
4930
4931     nnbl            = nbl_list->nnbl;
4932     nbl             = nbl_list->nbl;
4933     CombineNBLists  = nbl_list->bCombined;
4934
4935     if (debug)
4936     {
4937         fprintf(debug, "ns making %d nblists\n", nnbl);
4938     }
4939
4940     nbat->bUseBufferFlags = (nbat->nout > 1);
4941     /* We should re-init the flags before making the first list */
4942     if (nbat->bUseBufferFlags && (LOCAL_I(iloc) || bGPUCPU))
4943     {
4944         init_buffer_flags(&nbat->buffer_flags, nbat->natoms);
4945     }
4946
4947     if (nbl_list->bSimple)
4948     {
4949         switch (nb_kernel_type)
4950         {
4951 #ifdef GMX_NBNXN_SIMD_4XN
4952             case nbnxnk4xN_SIMD_4xN:
4953                 nbs->icell_set_x = icell_set_x_simd_4xn;
4954                 break;
4955 #endif
4956 #ifdef GMX_NBNXN_SIMD_2XNN
4957             case nbnxnk4xN_SIMD_2xNN:
4958                 nbs->icell_set_x = icell_set_x_simd_2xnn;
4959                 break;
4960 #endif
4961             default:
4962                 nbs->icell_set_x = icell_set_x_simple;
4963                 break;
4964         }
4965     }
4966     else
4967     {
4968 #ifdef NBNXN_SEARCH_BB_SIMD4
4969         nbs->icell_set_x = icell_set_x_supersub_simd4;
4970 #else
4971         nbs->icell_set_x = icell_set_x_supersub;
4972 #endif
4973     }
4974
4975     if (LOCAL_I(iloc))
4976     {
4977         /* Only zone (grid) 0 vs 0 */
4978         nzi = 1;
4979         zj0 = 0;
4980         zj1 = 1;
4981     }
4982     else
4983     {
4984         nzi = nbs->zones->nizone;
4985     }
4986
4987     if (!nbl_list->bSimple && min_ci_balanced > 0)
4988     {
4989         nsubpair_max = get_nsubpair_max(nbs, iloc, rlist, min_ci_balanced);
4990     }
4991     else
4992     {
4993         nsubpair_max = 0;
4994     }
4995
4996     /* Clear all pair-lists */
4997     for (th = 0; th < nnbl; th++)
4998     {
4999         clear_pairlist(nbl[th]);
5000     }
5001
5002     for (zi = 0; zi < nzi; zi++)
5003     {
5004         gridi = &nbs->grid[zi];
5005
5006         if (NONLOCAL_I(iloc))
5007         {
5008             zj0 = nbs->zones->izone[zi].j0;
5009             zj1 = nbs->zones->izone[zi].j1;
5010             if (zi == 0)
5011             {
5012                 zj0++;
5013             }
5014         }
5015         for (zj = zj0; zj < zj1; zj++)
5016         {
5017             gridj = &nbs->grid[zj];
5018
5019             if (debug)
5020             {
5021                 fprintf(debug, "ns search grid %d vs %d\n", zi, zj);
5022             }
5023
5024             nbs_cycle_start(&nbs->cc[enbsCCsearch]);
5025
5026             if (nbl[0]->bSimple && !gridi->bSimple)
5027             {
5028                 /* Hybrid list, determine blocking later */
5029                 ci_block = 0;
5030             }
5031             else
5032             {
5033                 ci_block = get_ci_block_size(gridi, nbs->DomDec, nnbl);
5034             }
5035
5036 #pragma omp parallel for num_threads(nnbl) schedule(static)
5037             for (th = 0; th < nnbl; th++)
5038             {
5039                 /* Re-init the thread-local work flag data before making
5040                  * the first list (not an elegant conditional).
5041                  */
5042                 if (nbat->bUseBufferFlags && ((zi == 0 && zj == 0) ||
5043                                               (bGPUCPU && zi == 0 && zj == 1)))
5044                 {
5045                     init_buffer_flags(&nbs->work[th].buffer_flags, nbat->natoms);
5046                 }
5047
5048                 if (CombineNBLists && th > 0)
5049                 {
5050                     clear_pairlist(nbl[th]);
5051                 }
5052
5053                 /* With GPU: generate progressively smaller lists for
5054                  * load balancing for local only or non-local with 2 zones.
5055                  */
5056                 progBal = (LOCAL_I(iloc) || nbs->zones->n <= 2);
5057
5058                 /* Divide the i super cell equally over the nblists */
5059                 nbnxn_make_pairlist_part(nbs, gridi, gridj,
5060                                          &nbs->work[th], nbat, excl,
5061                                          rlist,
5062                                          nb_kernel_type,
5063                                          ci_block,
5064                                          nbat->bUseBufferFlags,
5065                                          nsubpair_max,
5066                                          progBal, min_ci_balanced,
5067                                          th, nnbl,
5068                                          nbl[th]);
5069             }
5070             nbs_cycle_stop(&nbs->cc[enbsCCsearch]);
5071
5072             np_tot = 0;
5073             np_noq = 0;
5074             np_hlj = 0;
5075             for (th = 0; th < nnbl; th++)
5076             {
5077                 inc_nrnb(nrnb, eNR_NBNXN_DIST2, nbs->work[th].ndistc);
5078
5079                 if (nbl_list->bSimple)
5080                 {
5081                     np_tot += nbl[th]->ncj;
5082                     np_noq += nbl[th]->work->ncj_noq;
5083                     np_hlj += nbl[th]->work->ncj_hlj;
5084                 }
5085                 else
5086                 {
5087                     /* This count ignores potential subsequent pair pruning */
5088                     np_tot += nbl[th]->nci_tot;
5089                 }
5090             }
5091             nap                   = nbl[0]->na_ci*nbl[0]->na_cj;
5092             nbl_list->natpair_ljq = (np_tot - np_noq)*nap - np_hlj*nap/2;
5093             nbl_list->natpair_lj  = np_noq*nap;
5094             nbl_list->natpair_q   = np_hlj*nap/2;
5095
5096             if (CombineNBLists && nnbl > 1)
5097             {
5098                 nbs_cycle_start(&nbs->cc[enbsCCcombine]);
5099
5100                 combine_nblists(nnbl-1, nbl+1, nbl[0]);
5101
5102                 nbs_cycle_stop(&nbs->cc[enbsCCcombine]);
5103             }
5104         }
5105     }
5106
5107     if (!nbl_list->bSimple)
5108     {
5109         /* Sort the entries on size, large ones first */
5110         if (CombineNBLists || nnbl == 1)
5111         {
5112             sort_sci(nbl[0]);
5113         }
5114         else
5115         {
5116 #pragma omp parallel for num_threads(nnbl) schedule(static)
5117             for (th = 0; th < nnbl; th++)
5118             {
5119                 sort_sci(nbl[th]);
5120             }
5121         }
5122     }
5123
5124     if (nbat->bUseBufferFlags)
5125     {
5126         reduce_buffer_flags(nbs, nnbl, &nbat->buffer_flags);
5127     }
5128
5129     /* Special performance logging stuff (env.var. GMX_NBNXN_CYCLE) */
5130     if (LOCAL_I(iloc))
5131     {
5132         nbs->search_count++;
5133     }
5134     if (nbs->print_cycles &&
5135         (!nbs->DomDec || (nbs->DomDec && !LOCAL_I(iloc))) &&
5136         nbs->search_count % 100 == 0)
5137     {
5138         nbs_cycle_print(stderr, nbs);
5139     }
5140
5141     if (debug && (CombineNBLists && nnbl > 1))
5142     {
5143         if (nbl[0]->bSimple)
5144         {
5145             print_nblist_statistics_simple(debug, nbl[0], nbs, rlist);
5146         }
5147         else
5148         {
5149             print_nblist_statistics_supersub(debug, nbl[0], nbs, rlist);
5150         }
5151     }
5152
5153     if (debug)
5154     {
5155         if (gmx_debug_at)
5156         {
5157             if (nbl[0]->bSimple)
5158             {
5159                 print_nblist_ci_cj(debug, nbl[0]);
5160             }
5161             else
5162             {
5163                 print_nblist_sci_cj(debug, nbl[0]);
5164             }
5165         }
5166
5167         if (nbat->bUseBufferFlags)
5168         {
5169             print_reduction_cost(&nbat->buffer_flags, nnbl);
5170         }
5171     }
5172 }