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