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