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