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