src/gromacs/mdlib/nbnxn_atomdata.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 #include "gmxpre.h"
  37
  38 #include "nbnxn_atomdata.h"
  39
  40 #include <assert.h>
  41 #include <stdlib.h>
  42 #include <string.h>
  43
  44 #include <cmath>
  45
  46 #include <algorithm>
  47
  48 #include "thread_mpi/atomic.h"
  49
  50 #include "gromacs/math/functions.h"
  51 #include "gromacs/math/utilities.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  54 #include "gromacs/mdlib/nb_verlet.h"
  55 #include "gromacs/mdlib/nbnxn_consts.h"
  56 #include "gromacs/mdlib/nbnxn_internal.h"
  57 #include "gromacs/mdlib/nbnxn_search.h"
  58 #include "gromacs/mdlib/nbnxn_util.h"
  59 #include "gromacs/pbcutil/ishift.h"
  60 #include "gromacs/simd/simd.h"
  61 #include "gromacs/utility/exceptions.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/gmxomp.h"
  64 #include "gromacs/utility/smalloc.h"
  65
  66 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  67
  68 /* Default nbnxn allocation routine, allocates NBNXN_MEM_ALIGN byte aligned */
  69 void nbnxn_alloc_aligned(void **ptr, size_t nbytes)
  70 {
  71     *ptr = save_malloc_aligned("ptr", __FILE__, __LINE__, nbytes, 1, NBNXN_MEM_ALIGN);
  72 }
  73
  74 /* Free function for memory allocated with nbnxn_alloc_aligned */
  75 void nbnxn_free_aligned(void *ptr)
  76 {
  77     sfree_aligned(ptr);
  78 }
  79
  80 /* Reallocation wrapper function for nbnxn data structures */
  81 void nbnxn_realloc_void(void **ptr,
  82                         int nbytes_copy, int nbytes_new,
  83                         nbnxn_alloc_t *ma,
  84                         nbnxn_free_t  *mf)
  85 {
  86     void *ptr_new;
  87
  88     ma(&ptr_new, nbytes_new);
  89
  90     if (nbytes_new > 0 && ptr_new == nullptr)
  91     {
  92         gmx_fatal(FARGS, "Allocation of %d bytes failed", nbytes_new);
  93     }
  94
  95     if (nbytes_copy > 0)
  96     {
  97         if (nbytes_new < nbytes_copy)
  98         {
  99             gmx_incons("In nbnxn_realloc_void: new size less than copy size");
 100         }
 101         memcpy(ptr_new, *ptr, nbytes_copy);
 102     }
 103     if (*ptr != nullptr)
 104     {
 105         mf(*ptr);
 106     }
 107     *ptr = ptr_new;
 108 }
 109
 110 /* Reallocate the nbnxn_atomdata_t for a size of n atoms */
 111 void nbnxn_atomdata_realloc(nbnxn_atomdata_t *nbat, int n)
 112 {
 113     int t;
 114
 115     nbnxn_realloc_void((void **)&nbat->type,
 116                        nbat->natoms*sizeof(*nbat->type),
 117                        n*sizeof(*nbat->type),
 118                        nbat->alloc, nbat->free);
 119     nbnxn_realloc_void((void **)&nbat->lj_comb,
 120                        nbat->natoms*2*sizeof(*nbat->lj_comb),
 121                        n*2*sizeof(*nbat->lj_comb),
 122                        nbat->alloc, nbat->free);
 123     if (nbat->XFormat != nbatXYZQ)
 124     {
 125         nbnxn_realloc_void((void **)&nbat->q,
 126                            nbat->natoms*sizeof(*nbat->q),
 127                            n*sizeof(*nbat->q),
 128                            nbat->alloc, nbat->free);
 129     }
 130     if (nbat->nenergrp > 1)
 131     {
 132         nbnxn_realloc_void((void **)&nbat->energrp,
 133                            nbat->natoms/nbat->na_c*sizeof(*nbat->energrp),
 134                            n/nbat->na_c*sizeof(*nbat->energrp),
 135                            nbat->alloc, nbat->free);
 136     }
 137     nbnxn_realloc_void((void **)&nbat->x,
 138                        nbat->natoms*nbat->xstride*sizeof(*nbat->x),
 139                        n*nbat->xstride*sizeof(*nbat->x),
 140                        nbat->alloc, nbat->free);
 141     for (t = 0; t < nbat->nout; t++)
 142     {
 143         /* Allocate one element extra for possible signaling with GPUs */
 144         nbnxn_realloc_void((void **)&nbat->out[t].f,
 145                            nbat->natoms*nbat->fstride*sizeof(*nbat->out[t].f),
 146                            n*nbat->fstride*sizeof(*nbat->out[t].f),
 147                            nbat->alloc, nbat->free);
 148     }
 149     nbat->nalloc = n;
 150 }
 151
 152 /* Initializes an nbnxn_atomdata_output_t data structure */
 153 static void nbnxn_atomdata_output_init(nbnxn_atomdata_output_t *out,
 154                                        int nb_kernel_type,
 155                                        int nenergrp, int stride,
 156                                        nbnxn_alloc_t *ma)
 157 {
 158     out->f = nullptr;
 159     ma((void **)&out->fshift, SHIFTS*DIM*sizeof(*out->fshift));
 160     out->nV = nenergrp*nenergrp;
 161     ma((void **)&out->Vvdw, out->nV*sizeof(*out->Vvdw));
 162     ma((void **)&out->Vc, out->nV*sizeof(*out->Vc  ));
 163
 164     if (nb_kernel_type == nbnxnk4xN_SIMD_4xN ||
 165         nb_kernel_type == nbnxnk4xN_SIMD_2xNN)
 166     {
 167         int cj_size  = nbnxn_kernel_to_cluster_j_size(nb_kernel_type);
 168         out->nVS = nenergrp*nenergrp*stride*(cj_size>>1)*cj_size;
 169         ma((void **)&out->VSvdw, out->nVS*sizeof(*out->VSvdw));
 170         ma((void **)&out->VSc, out->nVS*sizeof(*out->VSc  ));
 171     }
 172     else
 173     {
 174         out->nVS = 0;
 175     }
 176 }
 177
 178 static void copy_int_to_nbat_int(const int *a, int na, int na_round,
 179                                  const int *in, int fill, int *innb)
 180 {
 181     int i, j;
 182
 183     j = 0;
 184     for (i = 0; i < na; i++)
 185     {
 186         innb[j++] = in[a[i]];
 187     }
 188     /* Complete the partially filled last cell with fill */
 189     for (; i < na_round; i++)
 190     {
 191         innb[j++] = fill;
 192     }
 193 }
 194
 195 void copy_rvec_to_nbat_real(const int *a, int na, int na_round,
 196                             const rvec *x, int nbatFormat,
 197                             real *xnb, int a0)
 198 {
 199     /* We complete partially filled cells, can only be the last one in each
 200      * column, with coordinates farAway. The actual coordinate value does
 201      * not influence the results, since these filler particles do not interact.
 202      * Clusters with normal atoms + fillers have a bounding box based only
 203      * on the coordinates of the atoms. Clusters with only fillers have as
 204      * the bounding box the coordinates of the first filler. Such clusters
 205      * are not considered as i-entries, but they are considered as j-entries.
 206      * So for performance it is better to have their bounding boxes far away,
 207      * such that filler only clusters don't end up in the pair list.
 208      */
 209     const real farAway = -1000000;
 210
 211     int        i, j, c;
 212
 213     switch (nbatFormat)
 214     {
 215         case nbatXYZ:
 216             j = a0*STRIDE_XYZ;
 217             for (i = 0; i < na; i++)
 218             {
 219                 xnb[j++] = x[a[i]][XX];
 220                 xnb[j++] = x[a[i]][YY];
 221                 xnb[j++] = x[a[i]][ZZ];
 222             }
 223             /* Complete the partially filled last cell with farAway elements */
 224             for (; i < na_round; i++)
 225             {
 226                 xnb[j++] = farAway;
 227                 xnb[j++] = farAway;
 228                 xnb[j++] = farAway;
 229             }
 230             break;
 231         case nbatXYZQ:
 232             j = a0*STRIDE_XYZQ;
 233             for (i = 0; i < na; i++)
 234             {
 235                 xnb[j++] = x[a[i]][XX];
 236                 xnb[j++] = x[a[i]][YY];
 237                 xnb[j++] = x[a[i]][ZZ];
 238                 j++;
 239             }
 240             /* Complete the partially filled last cell with zeros */
 241             for (; i < na_round; i++)
 242             {
 243                 xnb[j++] = farAway;
 244                 xnb[j++] = farAway;
 245                 xnb[j++] = farAway;
 246                 j++;
 247             }
 248             break;
 249         case nbatX4:
 250             j = atom_to_x_index<c_packX4>(a0);
 251             c = a0 & (c_packX4-1);
 252             for (i = 0; i < na; i++)
 253             {
 254                 xnb[j+XX*c_packX4] = x[a[i]][XX];
 255                 xnb[j+YY*c_packX4] = x[a[i]][YY];
 256                 xnb[j+ZZ*c_packX4] = x[a[i]][ZZ];
 257                 j++;
 258                 c++;
 259                 if (c == c_packX4)
 260                 {
 261                     j += (DIM-1)*c_packX4;
 262                     c  = 0;
 263                 }
 264             }
 265             /* Complete the partially filled last cell with zeros */
 266             for (; i < na_round; i++)
 267             {
 268                 xnb[j+XX*c_packX4] = farAway;
 269                 xnb[j+YY*c_packX4] = farAway;
 270                 xnb[j+ZZ*c_packX4] = farAway;
 271                 j++;
 272                 c++;
 273                 if (c == c_packX4)
 274                 {
 275                     j += (DIM-1)*c_packX4;
 276                     c  = 0;
 277                 }
 278             }
 279             break;
 280         case nbatX8:
 281             j = atom_to_x_index<c_packX8>(a0);
 282             c = a0 & (c_packX8 - 1);
 283             for (i = 0; i < na; i++)
 284             {
 285                 xnb[j+XX*c_packX8] = x[a[i]][XX];
 286                 xnb[j+YY*c_packX8] = x[a[i]][YY];
 287                 xnb[j+ZZ*c_packX8] = x[a[i]][ZZ];
 288                 j++;
 289                 c++;
 290                 if (c == c_packX8)
 291                 {
 292                     j += (DIM-1)*c_packX8;
 293                     c  = 0;
 294                 }
 295             }
 296             /* Complete the partially filled last cell with zeros */
 297             for (; i < na_round; i++)
 298             {
 299                 xnb[j+XX*c_packX8] = farAway;
 300                 xnb[j+YY*c_packX8] = farAway;
 301                 xnb[j+ZZ*c_packX8] = farAway;
 302                 j++;
 303                 c++;
 304                 if (c == c_packX8)
 305                 {
 306                     j += (DIM-1)*c_packX8;
 307                     c  = 0;
 308                 }
 309             }
 310             break;
 311         default:
 312             gmx_incons("Unsupported nbnxn_atomdata_t format");
 313     }
 314 }
 315
 316 /* Stores the LJ parameter data in a format convenient for different kernels */
 317 static void set_lj_parameter_data(nbnxn_atomdata_t *nbat, gmx_bool bSIMD)
 318 {
 319     real c6, c12;
 320
 321     int  nt = nbat->ntype;
 322
 323     if (bSIMD)
 324     {
 325 #if GMX_SIMD
 326         /* nbfp_aligned stores two parameters using the stride most suitable
 327          * for the present SIMD architecture, as specified by the constant
 328          * c_simdBestPairAlignment from the SIMD header.
 329          * There's a slight inefficiency in allocating and initializing nbfp_aligned
 330          * when it might not be used, but introducing the conditional code is not
 331          * really worth it.
 332          */
 333         nbat->alloc((void **)&nbat->nbfp_aligned,
 334                     nt*nt*c_simdBestPairAlignment*sizeof(*nbat->nbfp_aligned));
 335         for (int i = 0; i < nt; i++)
 336         {
 337             for (int j = 0; j < nt; j++)
 338             {
 339                 nbat->nbfp_aligned[(i*nt+j)*c_simdBestPairAlignment+0] = nbat->nbfp[(i*nt+j)*2+0];
 340                 nbat->nbfp_aligned[(i*nt+j)*c_simdBestPairAlignment+1] = nbat->nbfp[(i*nt+j)*2+1];
 341                 nbat->nbfp_aligned[(i*nt+j)*c_simdBestPairAlignment+2] = 0;
 342                 nbat->nbfp_aligned[(i*nt+j)*c_simdBestPairAlignment+3] = 0;
 343             }
 344         }
 345 #endif
 346     }
 347
 348     /* We use combination rule data for SIMD combination rule kernels
 349      * and with LJ-PME kernels. We then only need parameters per atom type,
 350      * not per pair of atom types.
 351      */
 352     switch (nbat->comb_rule)
 353     {
 354         case ljcrGEOM:
 355             nbat->comb_rule = ljcrGEOM;
 356
 357             for (int i = 0; i < nt; i++)
 358             {
 359                 /* Store the sqrt of the diagonal from the nbfp matrix */
 360                 nbat->nbfp_comb[i*2  ] = std::sqrt(nbat->nbfp[(i*nt+i)*2  ]);
 361                 nbat->nbfp_comb[i*2+1] = std::sqrt(nbat->nbfp[(i*nt+i)*2+1]);
 362             }
 363             break;
 364         case ljcrLB:
 365             for (int i = 0; i < nt; i++)
 366             {
 367                 /* Get 6*C6 and 12*C12 from the diagonal of the nbfp matrix */
 368                 c6  = nbat->nbfp[(i*nt+i)*2  ];
 369                 c12 = nbat->nbfp[(i*nt+i)*2+1];
 370                 if (c6 > 0 && c12 > 0)
 371                 {
 372                     /* We store 0.5*2^1/6*sigma and sqrt(4*3*eps),
 373                      * so we get 6*C6 and 12*C12 after combining.
 374                      */
 375                     nbat->nbfp_comb[i*2  ] = 0.5*gmx::sixthroot(c12/c6);
 376                     nbat->nbfp_comb[i*2+1] = std::sqrt(c6*c6/c12);
 377                 }
 378                 else
 379                 {
 380                     nbat->nbfp_comb[i*2  ] = 0;
 381                     nbat->nbfp_comb[i*2+1] = 0;
 382                 }
 383             }
 384             break;
 385         case ljcrNONE:
 386             /* We always store the full matrix (see code above) */
 387             break;
 388         default:
 389             gmx_incons("Unknown combination rule");
 390             break;
 391     }
 392 }
 393
 394 #if GMX_SIMD
 395 static void
 396 nbnxn_atomdata_init_simple_exclusion_masks(nbnxn_atomdata_t *nbat)
 397 {
 398     const int simd_width = GMX_SIMD_REAL_WIDTH;
 399     int       simd_excl_size;
 400     /* Set the diagonal cluster pair exclusion mask setup data.
 401      * In the kernel we check 0 < j - i to generate the masks.
 402      * Here we store j - i for generating the mask for the first i,
 403      * we substract 0.5 to avoid rounding issues.
 404      * In the kernel we can subtract 1 to generate the subsequent mask.
 405      */
 406     int        simd_4xn_diag_size;
 407
 408     simd_4xn_diag_size = std::max(NBNXN_CPU_CLUSTER_I_SIZE, simd_width);
 409     snew_aligned(nbat->simd_4xn_diagonal_j_minus_i, simd_4xn_diag_size, NBNXN_MEM_ALIGN);
 410     for (int j = 0; j < simd_4xn_diag_size; j++)
 411     {
 412         nbat->simd_4xn_diagonal_j_minus_i[j] = j - 0.5;
 413     }
 414
 415     snew_aligned(nbat->simd_2xnn_diagonal_j_minus_i, simd_width, NBNXN_MEM_ALIGN);
 416     for (int j = 0; j < simd_width/2; j++)
 417     {
 418         /* The j-cluster size is half the SIMD width */
 419         nbat->simd_2xnn_diagonal_j_minus_i[j]              = j - 0.5;
 420         /* The next half of the SIMD width is for i + 1 */
 421         nbat->simd_2xnn_diagonal_j_minus_i[simd_width/2+j] = j - 1 - 0.5;
 422     }
 423
 424     /* We use up to 32 bits for exclusion masking.
 425      * The same masks are used for the 4xN and 2x(N+N) kernels.
 426      * The masks are read either into integer SIMD registers or into
 427      * real SIMD registers (together with a cast).
 428      * In single precision this means the real and integer SIMD registers
 429      * are of equal size.
 430      */
 431     simd_excl_size = NBNXN_CPU_CLUSTER_I_SIZE*simd_width;
 432 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
 433     snew_aligned(nbat->simd_exclusion_filter64, simd_excl_size,   NBNXN_MEM_ALIGN);
 434 #else
 435     snew_aligned(nbat->simd_exclusion_filter, simd_excl_size,   NBNXN_MEM_ALIGN);
 436 #endif
 437
 438     for (int j = 0; j < simd_excl_size; j++)
 439     {
 440         /* Set the consecutive bits for masking pair exclusions */
 441 #if GMX_DOUBLE && !GMX_SIMD_HAVE_INT32_LOGICAL
 442         nbat->simd_exclusion_filter64[j]     = (1U << j);
 443 #else
 444         nbat->simd_exclusion_filter[j]       = (1U << j);
 445 #endif
 446     }
 447
 448 #if !GMX_SIMD_HAVE_LOGICAL && !GMX_SIMD_HAVE_INT32_LOGICAL
 449     // If the SIMD implementation has no bitwise logical operation support
 450     // whatsoever we cannot use the normal masking. Instead,
 451     // we generate a vector of all 2^4 possible ways an i atom
 452     // interacts with its 4 j atoms. Each array entry contains
 453     // GMX_SIMD_REAL_WIDTH values that are read with a single aligned SIMD load.
 454     // Since there is no logical value representation in this case, we use
 455     // any nonzero value to indicate 'true', while zero mean 'false'.
 456     // This can then be converted to a SIMD boolean internally in the SIMD
 457     // module by comparing to zero.
 458     // Each array entry encodes how this i atom will interact with the 4 j atoms.
 459     // Matching code exists in set_ci_top_excls() to generate indices into this array.
 460     // Those indices are used in the kernels.
 461
 462     simd_excl_size = NBNXN_CPU_CLUSTER_I_SIZE*NBNXN_CPU_CLUSTER_I_SIZE;
 463     const real simdFalse =  0.0;
 464     const real simdTrue  =  1.0;
 465     real      *simd_interaction_array;
 466
 467     snew_aligned(simd_interaction_array, simd_excl_size * GMX_SIMD_REAL_WIDTH, NBNXN_MEM_ALIGN);
 468     for (int j = 0; j < simd_excl_size; j++)
 469     {
 470         int index = j * GMX_SIMD_REAL_WIDTH;
 471         for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 472         {
 473             simd_interaction_array[index + i] = (j & (1 << i)) ? simdTrue : simdFalse;
 474         }
 475     }
 476     nbat->simd_interaction_array = simd_interaction_array;
 477 #endif
 478 }
 479 #endif
 480
 481 /* Initializes an nbnxn_atomdata_t data structure */
 482 void nbnxn_atomdata_init(FILE *fp,
 483                          nbnxn_atomdata_t *nbat,
 484                          int nb_kernel_type,
 485                          int enbnxninitcombrule,
 486                          int ntype, const real *nbfp,
 487                          int n_energygroups,
 488                          int nout,
 489                          nbnxn_alloc_t *alloc,
 490                          nbnxn_free_t  *free)
 491 {
 492     int      nth;
 493     real     c6, c12, tol;
 494     char    *ptr;
 495     gmx_bool simple, bCombGeom, bCombLB, bSIMD;
 496
 497     if (alloc == nullptr)
 498     {
 499         nbat->alloc = nbnxn_alloc_aligned;
 500     }
 501     else
 502     {
 503         nbat->alloc = alloc;
 504     }
 505     if (free == nullptr)
 506     {
 507         nbat->free = nbnxn_free_aligned;
 508     }
 509     else
 510     {
 511         nbat->free = free;
 512     }
 513
 514     if (debug)
 515     {
 516         fprintf(debug, "There are %d atom types in the system, adding one for nbnxn_atomdata_t\n", ntype);
 517     }
 518     nbat->ntype = ntype + 1;
 519     nbat->alloc((void **)&nbat->nbfp,
 520                 nbat->ntype*nbat->ntype*2*sizeof(*nbat->nbfp));
 521     nbat->alloc((void **)&nbat->nbfp_comb, nbat->ntype*2*sizeof(*nbat->nbfp_comb));
 522
 523     /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
 524      * force-field floating point parameters.
 525      */
 526     tol = 1e-5;
 527     ptr = getenv("GMX_LJCOMB_TOL");
 528     if (ptr != nullptr)
 529     {
 530         double dbl;
 531
 532         sscanf(ptr, "%lf", &dbl);
 533         tol = dbl;
 534     }
 535     bCombGeom = TRUE;
 536     bCombLB   = TRUE;
 537
 538     /* Temporarily fill nbat->nbfp_comb with sigma and epsilon
 539      * to check for the LB rule.
 540      */
 541     for (int i = 0; i < ntype; i++)
 542     {
 543         c6  = nbfp[(i*ntype+i)*2  ]/6.0;
 544         c12 = nbfp[(i*ntype+i)*2+1]/12.0;
 545         if (c6 > 0 && c12 > 0)
 546         {
 547             nbat->nbfp_comb[i*2  ] = gmx::sixthroot(c12/c6);
 548             nbat->nbfp_comb[i*2+1] = 0.25*c6*c6/c12;
 549         }
 550         else if (c6 == 0 && c12 == 0)
 551         {
 552             nbat->nbfp_comb[i*2  ] = 0;
 553             nbat->nbfp_comb[i*2+1] = 0;
 554         }
 555         else
 556         {
 557             /* Can not use LB rule with only dispersion or repulsion */
 558             bCombLB = FALSE;
 559         }
 560     }
 561
 562     for (int i = 0; i < nbat->ntype; i++)
 563     {
 564         for (int j = 0; j < nbat->ntype; j++)
 565         {
 566             if (i < ntype && j < ntype)
 567             {
 568                 /* fr->nbfp has been updated, so that array too now stores c6/c12 including
 569                  * the 6.0/12.0 prefactors to save 2 flops in the most common case (force-only).
 570                  */
 571                 c6  = nbfp[(i*ntype+j)*2  ];
 572                 c12 = nbfp[(i*ntype+j)*2+1];
 573                 nbat->nbfp[(i*nbat->ntype+j)*2  ] = c6;
 574                 nbat->nbfp[(i*nbat->ntype+j)*2+1] = c12;
 575
 576                 /* Compare 6*C6 and 12*C12 for geometric cobination rule */
 577                 bCombGeom = bCombGeom &&
 578                     gmx_within_tol(c6*c6, nbfp[(i*ntype+i)*2  ]*nbfp[(j*ntype+j)*2  ], tol) &&
 579                     gmx_within_tol(c12*c12, nbfp[(i*ntype+i)*2+1]*nbfp[(j*ntype+j)*2+1], tol);
 580
 581                 /* Compare C6 and C12 for Lorentz-Berthelot combination rule */
 582                 c6     /= 6.0;
 583                 c12    /= 12.0;
 584                 bCombLB = bCombLB &&
 585                     ((c6 == 0 && c12 == 0 &&
 586                       (nbat->nbfp_comb[i*2+1] == 0 || nbat->nbfp_comb[j*2+1] == 0)) ||
 587                      (c6 > 0 && c12 > 0 &&
 588                       gmx_within_tol(gmx::sixthroot(c12/c6),
 589                                      0.5*(nbat->nbfp_comb[i*2]+nbat->nbfp_comb[j*2]), tol) &&
 590                       gmx_within_tol(0.25*c6*c6/c12, std::sqrt(nbat->nbfp_comb[i*2+1]*nbat->nbfp_comb[j*2+1]), tol)));
 591             }
 592             else
 593             {
 594                 /* Add zero parameters for the additional dummy atom type */
 595                 nbat->nbfp[(i*nbat->ntype+j)*2  ] = 0;
 596                 nbat->nbfp[(i*nbat->ntype+j)*2+1] = 0;
 597             }
 598         }
 599     }
 600     if (debug)
 601     {
 602         fprintf(debug, "Combination rules: geometric %d Lorentz-Berthelot %d\n",
 603                 bCombGeom, bCombLB);
 604     }
 605
 606     simple = nbnxn_kernel_pairlist_simple(nb_kernel_type);
 607
 608     switch (enbnxninitcombrule)
 609     {
 610         case enbnxninitcombruleDETECT:
 611             /* We prefer the geometic combination rule,
 612              * as that gives a slightly faster kernel than the LB rule.
 613              */
 614             if (bCombGeom)
 615             {
 616                 nbat->comb_rule = ljcrGEOM;
 617             }
 618             else if (bCombLB)
 619             {
 620                 nbat->comb_rule = ljcrLB;
 621             }
 622             else
 623             {
 624                 nbat->comb_rule = ljcrNONE;
 625
 626                 nbat->free(nbat->nbfp_comb);
 627             }
 628
 629             if (fp)
 630             {
 631                 if (nbat->comb_rule == ljcrNONE)
 632                 {
 633                     fprintf(fp, "Using full Lennard-Jones parameter combination matrix\n\n");
 634                 }
 635                 else
 636                 {
 637                     fprintf(fp, "Using %s Lennard-Jones combination rule\n\n",
 638                             nbat->comb_rule == ljcrGEOM ? "geometric" : "Lorentz-Berthelot");
 639                 }
 640             }
 641             break;
 642         case enbnxninitcombruleGEOM:
 643             nbat->comb_rule = ljcrGEOM;
 644             break;
 645         case enbnxninitcombruleLB:
 646             nbat->comb_rule = ljcrLB;
 647             break;
 648         case enbnxninitcombruleNONE:
 649             nbat->comb_rule = ljcrNONE;
 650
 651             nbat->free(nbat->nbfp_comb);
 652             break;
 653         default:
 654             gmx_incons("Unknown enbnxninitcombrule");
 655     }
 656
 657     bSIMD = (nb_kernel_type == nbnxnk4xN_SIMD_4xN ||
 658              nb_kernel_type == nbnxnk4xN_SIMD_2xNN);
 659
 660     set_lj_parameter_data(nbat, bSIMD);
 661
 662     nbat->natoms  = 0;
 663     nbat->type    = nullptr;
 664     nbat->lj_comb = nullptr;
 665     if (simple)
 666     {
 667         int pack_x;
 668
 669         if (bSIMD)
 670         {
 671             pack_x = std::max(NBNXN_CPU_CLUSTER_I_SIZE,
 672                               nbnxn_kernel_to_cluster_j_size(nb_kernel_type));
 673             switch (pack_x)
 674             {
 675                 case 4:
 676                     nbat->XFormat = nbatX4;
 677                     break;
 678                 case 8:
 679                     nbat->XFormat = nbatX8;
 680                     break;
 681                 default:
 682                     gmx_incons("Unsupported packing width");
 683             }
 684         }
 685         else
 686         {
 687             nbat->XFormat = nbatXYZ;
 688         }
 689
 690         nbat->FFormat = nbat->XFormat;
 691     }
 692     else
 693     {
 694         nbat->XFormat = nbatXYZQ;
 695         nbat->FFormat = nbatXYZ;
 696     }
 697     nbat->q        = nullptr;
 698     nbat->nenergrp = n_energygroups;
 699     if (!simple)
 700     {
 701         /* Energy groups not supported yet for super-sub lists */
 702         if (n_energygroups > 1 && fp != nullptr)
 703         {
 704             fprintf(fp, "\nNOTE: With GPUs, reporting energy group contributions is not supported\n\n");
 705         }
 706         nbat->nenergrp = 1;
 707     }
 708     /* Temporary storage goes as #grp^3*simd_width^2/2, so limit to 64 */
 709     if (nbat->nenergrp > 64)
 710     {
 711         gmx_fatal(FARGS, "With NxN kernels not more than 64 energy groups are supported\n");
 712     }
 713     nbat->neg_2log = 1;
 714     while (nbat->nenergrp > (1<<nbat->neg_2log))
 715     {
 716         nbat->neg_2log++;
 717     }
 718     nbat->energrp = nullptr;
 719     nbat->alloc((void **)&nbat->shift_vec, SHIFTS*sizeof(*nbat->shift_vec));
 720     nbat->xstride = (nbat->XFormat == nbatXYZQ ? STRIDE_XYZQ : DIM);
 721     nbat->fstride = (nbat->FFormat == nbatXYZQ ? STRIDE_XYZQ : DIM);
 722     nbat->x       = nullptr;
 723
 724 #if GMX_SIMD
 725     if (simple)
 726     {
 727         nbnxn_atomdata_init_simple_exclusion_masks(nbat);
 728     }
 729 #endif
 730
 731     /* Initialize the output data structures */
 732     nbat->nout    = nout;
 733     snew(nbat->out, nbat->nout);
 734     nbat->nalloc  = 0;
 735     for (int i = 0; i < nbat->nout; i++)
 736     {
 737         nbnxn_atomdata_output_init(&nbat->out[i],
 738                                    nb_kernel_type,
 739                                    nbat->nenergrp, 1<<nbat->neg_2log,
 740                                    nbat->alloc);
 741     }
 742     nbat->buffer_flags.flag        = nullptr;
 743     nbat->buffer_flags.flag_nalloc = 0;
 744
 745     nth = gmx_omp_nthreads_get(emntNonbonded);
 746
 747     ptr = getenv("GMX_USE_TREEREDUCE");
 748     if (ptr != nullptr)
 749     {
 750         nbat->bUseTreeReduce = strtol(ptr, nullptr, 10);
 751     }
 752 #if defined __MIC__
 753     else if (nth > 8) /*on the CPU we currently don't benefit even at 32*/
 754     {
 755         nbat->bUseTreeReduce = 1;
 756     }
 757 #endif
 758     else
 759     {
 760         nbat->bUseTreeReduce = 0;
 761     }
 762     if (nbat->bUseTreeReduce)
 763     {
 764         if (fp)
 765         {
 766             fprintf(fp, "Using tree force reduction\n\n");
 767         }
 768         snew(nbat->syncStep, nth);
 769     }
 770 }
 771
 772 template<int packSize>
 773 static void copy_lj_to_nbat_lj_comb(const real *ljparam_type,
 774                                     const int *type, int na,
 775                                     real *ljparam_at)
 776 {
 777     /* The LJ params follow the combination rule:
 778      * copy the params for the type array to the atom array.
 779      */
 780     for (int is = 0; is < na; is += packSize)
 781     {
 782         for (int k = 0; k < packSize; k++)
 783         {
 784             int i = is + k;
 785             ljparam_at[is*2            + k] = ljparam_type[type[i]*2    ];
 786             ljparam_at[is*2 + packSize + k] = ljparam_type[type[i]*2 + 1];
 787         }
 788     }
 789 }
 790
 791 /* Sets the atom type in nbnxn_atomdata_t */
 792 static void nbnxn_atomdata_set_atomtypes(nbnxn_atomdata_t    *nbat,
 793                                          const nbnxn_search_t nbs,
 794                                          const int           *type)
 795 {
 796     for (int g = 0; g < nbs->ngrid; g++)
 797     {
 798         const nbnxn_grid_t * grid = &nbs->grid[g];
 799
 800         /* Loop over all columns and copy and fill */
 801         for (int i = 0; i < grid->ncx*grid->ncy; i++)
 802         {
 803             int ncz = grid->cxy_ind[i+1] - grid->cxy_ind[i];
 804             int ash = (grid->cell0 + grid->cxy_ind[i])*grid->na_sc;
 805
 806             copy_int_to_nbat_int(nbs->a+ash, grid->cxy_na[i], ncz*grid->na_sc,
 807                                  type, nbat->ntype-1, nbat->type+ash);
 808         }
 809     }
 810 }
 811
 812 /* Sets the LJ combination rule parameters in nbnxn_atomdata_t */
 813 static void nbnxn_atomdata_set_ljcombparams(nbnxn_atomdata_t    *nbat,
 814                                             const nbnxn_search_t nbs)
 815 {
 816     if (nbat->comb_rule != ljcrNONE)
 817     {
 818         for (int g = 0; g < nbs->ngrid; g++)
 819         {
 820             const nbnxn_grid_t * grid = &nbs->grid[g];
 821
 822             /* Loop over all columns and copy and fill */
 823             for (int i = 0; i < grid->ncx*grid->ncy; i++)
 824             {
 825                 int ncz = grid->cxy_ind[i+1] - grid->cxy_ind[i];
 826                 int ash = (grid->cell0 + grid->cxy_ind[i])*grid->na_sc;
 827
 828                 if (nbat->XFormat == nbatX4)
 829                 {
 830                     copy_lj_to_nbat_lj_comb<c_packX4>(nbat->nbfp_comb,
 831                                                       nbat->type + ash,
 832                                                       ncz*grid->na_sc,
 833                                                       nbat->lj_comb + ash*2);
 834                 }
 835                 else if (nbat->XFormat == nbatX8)
 836                 {
 837                     copy_lj_to_nbat_lj_comb<c_packX8>(nbat->nbfp_comb,
 838                                                       nbat->type + ash,
 839                                                       ncz*grid->na_sc,
 840                                                       nbat->lj_comb + ash*2);
 841                 }
 842                 else if (nbat->XFormat == nbatXYZQ)
 843                 {
 844                     copy_lj_to_nbat_lj_comb<1>(nbat->nbfp_comb,
 845                                                nbat->type + ash,
 846                                                ncz*grid->na_sc,
 847                                                nbat->lj_comb + ash*2);
 848                 }
 849             }
 850         }
 851     }
 852 }
 853
 854 /* Sets the charges in nbnxn_atomdata_t *nbat */
 855 static void nbnxn_atomdata_set_charges(nbnxn_atomdata_t    *nbat,
 856                                        const nbnxn_search_t nbs,
 857                                        const real          *charge)
 858 {
 859     int                 i;
 860     real               *q;
 861
 862     for (int g = 0; g < nbs->ngrid; g++)
 863     {
 864         const nbnxn_grid_t * grid = &nbs->grid[g];
 865
 866         /* Loop over all columns and copy and fill */
 867         for (int cxy = 0; cxy < grid->ncx*grid->ncy; cxy++)
 868         {
 869             int ash      = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
 870             int na       = grid->cxy_na[cxy];
 871             int na_round = (grid->cxy_ind[cxy+1] - grid->cxy_ind[cxy])*grid->na_sc;
 872
 873             if (nbat->XFormat == nbatXYZQ)
 874             {
 875                 q = nbat->x + ash*STRIDE_XYZQ + ZZ + 1;
 876                 for (i = 0; i < na; i++)
 877                 {
 878                     *q = charge[nbs->a[ash+i]];
 879                     q += STRIDE_XYZQ;
 880                 }
 881                 /* Complete the partially filled last cell with zeros */
 882                 for (; i < na_round; i++)
 883                 {
 884                     *q = 0;
 885                     q += STRIDE_XYZQ;
 886                 }
 887             }
 888             else
 889             {
 890                 q = nbat->q + ash;
 891                 for (i = 0; i < na; i++)
 892                 {
 893                     *q = charge[nbs->a[ash+i]];
 894                     q++;
 895                 }
 896                 /* Complete the partially filled last cell with zeros */
 897                 for (; i < na_round; i++)
 898                 {
 899                     *q = 0;
 900                     q++;
 901                 }
 902             }
 903         }
 904     }
 905 }
 906
 907 /* Set the charges of perturbed atoms in nbnxn_atomdata_t to 0.
 908  * This is to automatically remove the RF/PME self term in the nbnxn kernels.
 909  * Part of the zero interactions are still calculated in the normal kernels.
 910  * All perturbed interactions are calculated in the free energy kernel,
 911  * using the original charge and LJ data, not nbnxn_atomdata_t.
 912  */
 913 static void nbnxn_atomdata_mask_fep(nbnxn_atomdata_t    *nbat,
 914                                     const nbnxn_search_t nbs)
 915 {
 916     real               *q;
 917     int                 stride_q, nsubc;
 918
 919     if (nbat->XFormat == nbatXYZQ)
 920     {
 921         q        = nbat->x + ZZ + 1;
 922         stride_q = STRIDE_XYZQ;
 923     }
 924     else
 925     {
 926         q        = nbat->q;
 927         stride_q = 1;
 928     }
 929
 930     for (int g = 0; g < nbs->ngrid; g++)
 931     {
 932         const nbnxn_grid_t * grid = &nbs->grid[g];
 933         if (grid->bSimple)
 934         {
 935             nsubc = 1;
 936         }
 937         else
 938         {
 939             nsubc = c_gpuNumClusterPerCell;
 940         }
 941
 942         int c_offset = grid->cell0*grid->na_sc;
 943
 944         /* Loop over all columns and copy and fill */
 945         for (int c = 0; c < grid->nc*nsubc; c++)
 946         {
 947             /* Does this cluster contain perturbed particles? */
 948             if (grid->fep[c] != 0)
 949             {
 950                 for (int i = 0; i < grid->na_c; i++)
 951                 {
 952                     /* Is this a perturbed particle? */
 953                     if (grid->fep[c] & (1 << i))
 954                     {
 955                         int ind = c_offset + c*grid->na_c + i;
 956                         /* Set atom type and charge to non-interacting */
 957                         nbat->type[ind] = nbat->ntype - 1;
 958                         q[ind*stride_q] = 0;
 959                     }
 960                 }
 961             }
 962         }
 963     }
 964 }
 965
 966 /* Copies the energy group indices to a reordered and packed array */
 967 static void copy_egp_to_nbat_egps(const int *a, int na, int na_round,
 968                                   int na_c, int bit_shift,
 969                                   const int *in, int *innb)
 970 {
 971     int i;
 972     int comb;
 973
 974     int j = 0;
 975     for (i = 0; i < na; i += na_c)
 976     {
 977         /* Store na_c energy group numbers into one int */
 978         comb = 0;
 979         for (int sa = 0; sa < na_c; sa++)
 980         {
 981             int at = a[i+sa];
 982             if (at >= 0)
 983             {
 984                 comb |= (GET_CGINFO_GID(in[at]) << (sa*bit_shift));
 985             }
 986         }
 987         innb[j++] = comb;
 988     }
 989     /* Complete the partially filled last cell with fill */
 990     for (; i < na_round; i += na_c)
 991     {
 992         innb[j++] = 0;
 993     }
 994 }
 995
 996 /* Set the energy group indices for atoms in nbnxn_atomdata_t */
 997 static void nbnxn_atomdata_set_energygroups(nbnxn_atomdata_t    *nbat,
 998                                             const nbnxn_search_t nbs,
 999                                             const int           *atinfo)
1000 {
1001     if (nbat->nenergrp == 1)
1002     {
1003         return;
1004     }
1005
1006     for (int g = 0; g < nbs->ngrid; g++)
1007     {
1008         const nbnxn_grid_t * grid = &nbs->grid[g];
1009
1010         /* Loop over all columns and copy and fill */
1011         for (int i = 0; i < grid->ncx*grid->ncy; i++)
1012         {
1013             int ncz = grid->cxy_ind[i+1] - grid->cxy_ind[i];
1014             int ash = (grid->cell0 + grid->cxy_ind[i])*grid->na_sc;
1015
1016             copy_egp_to_nbat_egps(nbs->a+ash, grid->cxy_na[i], ncz*grid->na_sc,
1017                                   nbat->na_c, nbat->neg_2log,
1018                                   atinfo, nbat->energrp+(ash>>grid->na_c_2log));
1019         }
1020     }
1021 }
1022
1023 /* Sets all required atom parameter data in nbnxn_atomdata_t */
1024 void nbnxn_atomdata_set(nbnxn_atomdata_t    *nbat,
1025                         const nbnxn_search_t nbs,
1026                         const t_mdatoms     *mdatoms,
1027                         const int           *atinfo)
1028 {
1029     nbnxn_atomdata_set_atomtypes(nbat, nbs, mdatoms->typeA);
1030
1031     nbnxn_atomdata_set_charges(nbat, nbs, mdatoms->chargeA);
1032
1033     if (nbs->bFEP)
1034     {
1035         nbnxn_atomdata_mask_fep(nbat, nbs);
1036     }
1037
1038     /* This must be done after masking types for FEP */
1039     nbnxn_atomdata_set_ljcombparams(nbat, nbs);
1040
1041     nbnxn_atomdata_set_energygroups(nbat, nbs, atinfo);
1042 }
1043
1044 /* Copies the shift vector array to nbnxn_atomdata_t */
1045 void nbnxn_atomdata_copy_shiftvec(gmx_bool          bDynamicBox,
1046                                   rvec             *shift_vec,
1047                                   nbnxn_atomdata_t *nbat)
1048 {
1049     int i;
1050
1051     nbat->bDynamicBox = bDynamicBox;
1052     for (i = 0; i < SHIFTS; i++)
1053     {
1054         copy_rvec(shift_vec[i], nbat->shift_vec[i]);
1055     }
1056 }
1057
1058 /* Copies (and reorders) the coordinates to nbnxn_atomdata_t */
1059 void nbnxn_atomdata_copy_x_to_nbat_x(const nbnxn_search_t nbs,
1060                                      int                  locality,
1061                                      gmx_bool             FillLocal,
1062                                      rvec                *x,
1063                                      nbnxn_atomdata_t    *nbat)
1064 {
1065     int g0 = 0, g1 = 0;
1066     int nth, th;
1067
1068     switch (locality)
1069     {
1070         case eatAll:
1071             g0 = 0;
1072             g1 = nbs->ngrid;
1073             break;
1074         case eatLocal:
1075             g0 = 0;
1076             g1 = 1;
1077             break;
1078         case eatNonlocal:
1079             g0 = 1;
1080             g1 = nbs->ngrid;
1081             break;
1082     }
1083
1084     if (FillLocal)
1085     {
1086         nbat->natoms_local = nbs->grid[0].nc*nbs->grid[0].na_sc;
1087     }
1088
1089     nth = gmx_omp_nthreads_get(emntPairsearch);
1090
1091 #pragma omp parallel for num_threads(nth) schedule(static)
1092     for (th = 0; th < nth; th++)
1093     {
1094         try
1095         {
1096             for (int g = g0; g < g1; g++)
1097             {
1098                 const nbnxn_grid_t *grid;
1099                 int                 cxy0, cxy1;
1100
1101                 grid = &nbs->grid[g];
1102
1103                 cxy0 = (grid->ncx*grid->ncy* th   +nth-1)/nth;
1104                 cxy1 = (grid->ncx*grid->ncy*(th+1)+nth-1)/nth;
1105
1106                 for (int cxy = cxy0; cxy < cxy1; cxy++)
1107                 {
1108                     int na, ash, na_fill;
1109
1110                     na  = grid->cxy_na[cxy];
1111                     ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
1112
1113                     if (g == 0 && FillLocal)
1114                     {
1115                         na_fill =
1116                             (grid->cxy_ind[cxy+1] - grid->cxy_ind[cxy])*grid->na_sc;
1117                     }
1118                     else
1119                     {
1120                         /* We fill only the real particle locations.
1121                          * We assume the filling entries at the end have been
1122                          * properly set before during pair-list generation.
1123                          */
1124                         na_fill = na;
1125                     }
1126                     copy_rvec_to_nbat_real(nbs->a+ash, na, na_fill, x,
1127                                            nbat->XFormat, nbat->x, ash);
1128                 }
1129             }
1130         }
1131         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1132     }
1133 }
1134
1135 static void
1136 nbnxn_atomdata_clear_reals(real * gmx_restrict dest,
1137                            int i0, int i1)
1138 {
1139     for (int i = i0; i < i1; i++)
1140     {
1141         dest[i] = 0;
1142     }
1143 }
1144
1145 gmx_unused static void
1146 nbnxn_atomdata_reduce_reals(real * gmx_restrict dest,
1147                             gmx_bool bDestSet,
1148                             real ** gmx_restrict src,
1149                             int nsrc,
1150                             int i0, int i1)
1151 {
1152     if (bDestSet)
1153     {
1154         /* The destination buffer contains data, add to it */
1155         for (int i = i0; i < i1; i++)
1156         {
1157             for (int s = 0; s < nsrc; s++)
1158             {
1159                 dest[i] += src[s][i];
1160             }
1161         }
1162     }
1163     else
1164     {
1165         /* The destination buffer is unitialized, set it first */
1166         for (int i = i0; i < i1; i++)
1167         {
1168             dest[i] = src[0][i];
1169             for (int s = 1; s < nsrc; s++)
1170             {
1171                 dest[i] += src[s][i];
1172             }
1173         }
1174     }
1175 }
1176
1177 gmx_unused static void
1178 nbnxn_atomdata_reduce_reals_simd(real gmx_unused * gmx_restrict dest,
1179                                  gmx_bool gmx_unused bDestSet,
1180                                  real gmx_unused ** gmx_restrict src,
1181                                  int gmx_unused nsrc,
1182                                  int gmx_unused i0, int gmx_unused i1)
1183 {
1184 #if GMX_SIMD
1185 /* The SIMD width here is actually independent of that in the kernels,
1186  * but we use the same width for simplicity (usually optimal anyhow).
1187  */
1188     SimdReal dest_SSE, src_SSE;
1189
1190     if (bDestSet)
1191     {
1192         for (int i = i0; i < i1; i += GMX_SIMD_REAL_WIDTH)
1193         {
1194             dest_SSE = load<SimdReal>(dest+i);
1195             for (int s = 0; s < nsrc; s++)
1196             {
1197                 src_SSE  = load<SimdReal>(src[s]+i);
1198                 dest_SSE = dest_SSE + src_SSE;
1199             }
1200             store(dest+i, dest_SSE);
1201         }
1202     }
1203     else
1204     {
1205         for (int i = i0; i < i1; i += GMX_SIMD_REAL_WIDTH)
1206         {
1207             dest_SSE = load<SimdReal>(src[0]+i);
1208             for (int s = 1; s < nsrc; s++)
1209             {
1210                 src_SSE  = load<SimdReal>(src[s]+i);
1211                 dest_SSE = dest_SSE + src_SSE;
1212             }
1213             store(dest+i, dest_SSE);
1214         }
1215     }
1216 #endif
1217 }
1218
1219 /* Add part of the force array(s) from nbnxn_atomdata_t to f */
1220 static void
1221 nbnxn_atomdata_add_nbat_f_to_f_part(const nbnxn_search_t nbs,
1222                                     const nbnxn_atomdata_t *nbat,
1223                                     nbnxn_atomdata_output_t *out,
1224                                     int nfa,
1225                                     int a0, int a1,
1226                                     rvec *f)
1227 {
1228     const int  *cell;
1229     const real *fnb;
1230
1231     cell = nbs->cell;
1232
1233     /* Loop over all columns and copy and fill */
1234     switch (nbat->FFormat)
1235     {
1236         case nbatXYZ:
1237         case nbatXYZQ:
1238             if (nfa == 1)
1239             {
1240                 fnb = out[0].f;
1241
1242                 for (int a = a0; a < a1; a++)
1243                 {
1244                     int i = cell[a]*nbat->fstride;
1245
1246                     f[a][XX] += fnb[i];
1247                     f[a][YY] += fnb[i+1];
1248                     f[a][ZZ] += fnb[i+2];
1249                 }
1250             }
1251             else
1252             {
1253                 for (int a = a0; a < a1; a++)
1254                 {
1255                     int i = cell[a]*nbat->fstride;
1256
1257                     for (int fa = 0; fa < nfa; fa++)
1258                     {
1259                         f[a][XX] += out[fa].f[i];
1260                         f[a][YY] += out[fa].f[i+1];
1261                         f[a][ZZ] += out[fa].f[i+2];
1262                     }
1263                 }
1264             }
1265             break;
1266         case nbatX4:
1267             if (nfa == 1)
1268             {
1269                 fnb = out[0].f;
1270
1271                 for (int a = a0; a < a1; a++)
1272                 {
1273                     int i = atom_to_x_index<c_packX4>(cell[a]);
1274
1275                     f[a][XX] += fnb[i+XX*c_packX4];
1276                     f[a][YY] += fnb[i+YY*c_packX4];
1277                     f[a][ZZ] += fnb[i+ZZ*c_packX4];
1278                 }
1279             }
1280             else
1281             {
1282                 for (int a = a0; a < a1; a++)
1283                 {
1284                     int i = atom_to_x_index<c_packX4>(cell[a]);
1285
1286                     for (int fa = 0; fa < nfa; fa++)
1287                     {
1288                         f[a][XX] += out[fa].f[i+XX*c_packX4];
1289                         f[a][YY] += out[fa].f[i+YY*c_packX4];
1290                         f[a][ZZ] += out[fa].f[i+ZZ*c_packX4];
1291                     }
1292                 }
1293             }
1294             break;
1295         case nbatX8:
1296             if (nfa == 1)
1297             {
1298                 fnb = out[0].f;
1299
1300                 for (int a = a0; a < a1; a++)
1301                 {
1302                     int i = atom_to_x_index<c_packX8>(cell[a]);
1303
1304                     f[a][XX] += fnb[i+XX*c_packX8];
1305                     f[a][YY] += fnb[i+YY*c_packX8];
1306                     f[a][ZZ] += fnb[i+ZZ*c_packX8];
1307                 }
1308             }
1309             else
1310             {
1311                 for (int a = a0; a < a1; a++)
1312                 {
1313                     int i = atom_to_x_index<c_packX8>(cell[a]);
1314
1315                     for (int fa = 0; fa < nfa; fa++)
1316                     {
1317                         f[a][XX] += out[fa].f[i+XX*c_packX8];
1318                         f[a][YY] += out[fa].f[i+YY*c_packX8];
1319                         f[a][ZZ] += out[fa].f[i+ZZ*c_packX8];
1320                     }
1321                 }
1322             }
1323             break;
1324         default:
1325             gmx_incons("Unsupported nbnxn_atomdata_t format");
1326     }
1327 }
1328
1329 static gmx_inline unsigned char reverse_bits(unsigned char b)
1330 {
1331     /* http://graphics.stanford.edu/~seander/bithacks.html#ReverseByteWith64BitsDiv */
1332     return (b * 0x0202020202ULL & 0x010884422010ULL) % 1023;
1333 }
1334
1335 static void nbnxn_atomdata_add_nbat_f_to_f_treereduce(const nbnxn_atomdata_t *nbat,
1336                                                       int                     nth)
1337 {
1338     const nbnxn_buffer_flags_t *flags = &nbat->buffer_flags;
1339
1340     int next_pow2 = 1<<(gmx::log2I(nth-1)+1);
1341
1342     assert(nbat->nout == nth); /* tree-reduce currently only works for nout==nth */
1343
1344     memset(nbat->syncStep, 0, sizeof(*(nbat->syncStep))*nth);
1345
1346 #pragma omp parallel num_threads(nth)
1347     {
1348         try
1349         {
1350             int   b0, b1, b;
1351             int   i0, i1;
1352             int   group_size, th;
1353
1354             th = gmx_omp_get_thread_num();
1355
1356             for (group_size = 2; group_size < 2*next_pow2; group_size *= 2)
1357             {
1358                 int index[2], group_pos, partner_pos, wu;
1359                 int partner_th = th ^ (group_size/2);
1360
1361                 if (group_size > 2)
1362                 {
1363 #ifdef TMPI_ATOMICS
1364                     /* wait on partner thread - replaces full barrier */
1365                     int sync_th, sync_group_size;
1366
1367                     tMPI_Atomic_memory_barrier();                         /* gurantee data is saved before marking work as done */
1368                     tMPI_Atomic_set(&(nbat->syncStep[th]), group_size/2); /* mark previous step as completed */
1369
1370                     /* find thread to sync with. Equal to partner_th unless nth is not a power of two. */
1371                     for (sync_th = partner_th, sync_group_size = group_size; sync_th >= nth && sync_group_size > 2; sync_group_size /= 2)
1372                     {
1373                         sync_th &= ~(sync_group_size/4);
1374                     }
1375                     if (sync_th < nth) /* otherwise nothing to sync index[1] will be >=nout */
1376                     {
1377                         /* wait on the thread which computed input data in previous step */
1378                         while (tMPI_Atomic_get((volatile tMPI_Atomic_t*)&(nbat->syncStep[sync_th])) < group_size/2)
1379                         {
1380                             gmx_pause();
1381                         }
1382                         /* guarantee that no later load happens before wait loop is finisehd */
1383                         tMPI_Atomic_memory_barrier();
1384                     }
1385 #else               /* TMPI_ATOMICS */
1386 #pragma omp barrier
1387 #endif
1388                 }
1389
1390                 /* Calculate buffers to sum (result goes into first buffer) */
1391                 group_pos = th % group_size;
1392                 index[0]  = th - group_pos;
1393                 index[1]  = index[0] + group_size/2;
1394
1395                 /* If no second buffer, nothing to do */
1396                 if (index[1] >= nbat->nout && group_size > 2)
1397                 {
1398                     continue;
1399                 }
1400
1401 #if NBNXN_BUFFERFLAG_MAX_THREADS > 256
1402 #error reverse_bits assumes max 256 threads
1403 #endif
1404                 /* Position is permuted so that one of the 2 vectors being added was computed on the same thread in the previous step.
1405                    This improves locality and enables to sync with just a single thread between steps (=the levels in the btree).
1406                    The permutation which allows this corresponds to reversing the bits of the group position.
1407                  */
1408                 group_pos = reverse_bits(group_pos)/(256/group_size);
1409
1410                 partner_pos = group_pos ^ 1;
1411
1412                 /* loop over two work-units (own and partner) */
1413                 for (wu = 0; wu < 2; wu++)
1414                 {
1415                     if (wu == 1)
1416                     {
1417                         if (partner_th < nth)
1418                         {
1419                             break; /* partner exists we don't have to do his work */
1420                         }
1421                         else
1422                         {
1423                             group_pos = partner_pos;
1424                         }
1425                     }
1426
1427                     /* Calculate the cell-block range for our thread */
1428                     b0 = (flags->nflag* group_pos   )/group_size;
1429                     b1 = (flags->nflag*(group_pos+1))/group_size;
1430
1431                     for (b = b0; b < b1; b++)
1432                     {
1433                         i0 =  b   *NBNXN_BUFFERFLAG_SIZE*nbat->fstride;
1434                         i1 = (b+1)*NBNXN_BUFFERFLAG_SIZE*nbat->fstride;
1435
1436                         if (bitmask_is_set(flags->flag[b], index[1]) || group_size > 2)
1437                         {
1438 #if GMX_SIMD
1439                             nbnxn_atomdata_reduce_reals_simd
1440 #else
1441                             nbnxn_atomdata_reduce_reals
1442 #endif
1443                                 (nbat->out[index[0]].f,
1444                                 bitmask_is_set(flags->flag[b], index[0]) || group_size > 2,
1445                                 &(nbat->out[index[1]].f), 1, i0, i1);
1446
1447                         }
1448                         else if (!bitmask_is_set(flags->flag[b], index[0]))
1449                         {
1450                             nbnxn_atomdata_clear_reals(nbat->out[index[0]].f,
1451                                                        i0, i1);
1452                         }
1453                     }
1454                 }
1455             }
1456         }
1457         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1458     }
1459 }
1460
1461
1462 static void nbnxn_atomdata_add_nbat_f_to_f_stdreduce(const nbnxn_atomdata_t *nbat,
1463                                                      int                     nth)
1464 {
1465 #pragma omp parallel for num_threads(nth) schedule(static)
1466     for (int th = 0; th < nth; th++)
1467     {
1468         try
1469         {
1470             const nbnxn_buffer_flags_t *flags;
1471             int   nfptr;
1472             real *fptr[NBNXN_BUFFERFLAG_MAX_THREADS];
1473
1474             flags = &nbat->buffer_flags;
1475
1476             /* Calculate the cell-block range for our thread */
1477             int b0 = (flags->nflag* th   )/nth;
1478             int b1 = (flags->nflag*(th+1))/nth;
1479
1480             for (int b = b0; b < b1; b++)
1481             {
1482                 int i0 =  b   *NBNXN_BUFFERFLAG_SIZE*nbat->fstride;
1483                 int i1 = (b+1)*NBNXN_BUFFERFLAG_SIZE*nbat->fstride;
1484
1485                 nfptr = 0;
1486                 for (int out = 1; out < nbat->nout; out++)
1487                 {
1488                     if (bitmask_is_set(flags->flag[b], out))
1489                     {
1490                         fptr[nfptr++] = nbat->out[out].f;
1491                     }
1492                 }
1493                 if (nfptr > 0)
1494                 {
1495 #if GMX_SIMD
1496                     nbnxn_atomdata_reduce_reals_simd
1497 #else
1498                     nbnxn_atomdata_reduce_reals
1499 #endif
1500                         (nbat->out[0].f,
1501                         bitmask_is_set(flags->flag[b], 0),
1502                         fptr, nfptr,
1503                         i0, i1);
1504                 }
1505                 else if (!bitmask_is_set(flags->flag[b], 0))
1506                 {
1507                     nbnxn_atomdata_clear_reals(nbat->out[0].f,
1508                                                i0, i1);
1509                 }
1510             }
1511         }
1512         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1513     }
1514 }
1515
1516 /* Add the force array(s) from nbnxn_atomdata_t to f */
1517 void nbnxn_atomdata_add_nbat_f_to_f(const nbnxn_search_t    nbs,
1518                                     int                     locality,
1519                                     const nbnxn_atomdata_t *nbat,
1520                                     rvec                   *f)
1521 {
1522     int a0 = 0, na = 0;
1523
1524     nbs_cycle_start(&nbs->cc[enbsCCreducef]);
1525
1526     switch (locality)
1527     {
1528         case eatAll:
1529             a0 = 0;
1530             na = nbs->natoms_nonlocal;
1531             break;
1532         case eatLocal:
1533             a0 = 0;
1534             na = nbs->natoms_local;
1535             break;
1536         case eatNonlocal:
1537             a0 = nbs->natoms_local;
1538             na = nbs->natoms_nonlocal - nbs->natoms_local;
1539             break;
1540     }
1541
1542     int nth = gmx_omp_nthreads_get(emntNonbonded);
1543
1544     if (nbat->nout > 1)
1545     {
1546         if (locality != eatAll)
1547         {
1548             gmx_incons("add_f_to_f called with nout>1 and locality!=eatAll");
1549         }
1550
1551         /* Reduce the force thread output buffers into buffer 0, before adding
1552          * them to the, differently ordered, "real" force buffer.
1553          */
1554         if (nbat->bUseTreeReduce)
1555         {
1556             nbnxn_atomdata_add_nbat_f_to_f_treereduce(nbat, nth);
1557         }
1558         else
1559         {
1560             nbnxn_atomdata_add_nbat_f_to_f_stdreduce(nbat, nth);
1561         }
1562     }
1563 #pragma omp parallel for num_threads(nth) schedule(static)
1564     for (int th = 0; th < nth; th++)
1565     {
1566         try
1567         {
1568             nbnxn_atomdata_add_nbat_f_to_f_part(nbs, nbat,
1569                                                 nbat->out,
1570                                                 1,
1571                                                 a0+((th+0)*na)/nth,
1572                                                 a0+((th+1)*na)/nth,
1573                                                 f);
1574         }
1575         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1576     }
1577
1578     nbs_cycle_stop(&nbs->cc[enbsCCreducef]);
1579 }
1580
1581 /* Adds the shift forces from nbnxn_atomdata_t to fshift */
1582 void nbnxn_atomdata_add_nbat_fshift_to_fshift(const nbnxn_atomdata_t *nbat,
1583                                               rvec                   *fshift)
1584 {
1585     const nbnxn_atomdata_output_t * out = nbat->out;
1586
1587     for (int s = 0; s < SHIFTS; s++)
1588     {
1589         rvec sum;
1590         clear_rvec(sum);
1591         for (int th = 0; th < nbat->nout; th++)
1592         {
1593             sum[XX] += out[th].fshift[s*DIM+XX];
1594             sum[YY] += out[th].fshift[s*DIM+YY];
1595             sum[ZZ] += out[th].fshift[s*DIM+ZZ];
1596         }
1597         rvec_inc(fshift[s], sum);
1598     }
1599 }