src/gromacs/mdlib/nbnxn_atomdata.c

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