src/mdlib/nbnxn_atomdata.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2012, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <math.h>
  44 #include <string.h>
  45 #include "smalloc.h"
  46 #include "macros.h"
  47 #include "vec.h"
  48 #include "nbnxn_consts.h"
  49 #include "nbnxn_internal.h"
  50 #include "nbnxn_search.h"
  51 #include "nbnxn_atomdata.h"
  52 #include "gmx_omp_nthreads.h"
  53
  54 /* Default nbnxn allocation routine, allocates 32 byte aligned,
  55  * which works for plain C and aligned SSE and AVX loads/stores.
  56  */
  57 void nbnxn_alloc_aligned(void **ptr,size_t nbytes)
  58 {
  59     *ptr = save_malloc_aligned("ptr",__FILE__,__LINE__,nbytes,1,32);
  60 }
  61
  62 /* Free function for memory allocated with nbnxn_alloc_aligned */
  63 void nbnxn_free_aligned(void *ptr)
  64 {
  65     sfree_aligned(ptr);
  66 }
  67
  68 /* Reallocation wrapper function for nbnxn data structures */
  69 void nbnxn_realloc_void(void **ptr,
  70                         int nbytes_copy,int nbytes_new,
  71                         nbnxn_alloc_t *ma,
  72                         nbnxn_free_t  *mf)
  73 {
  74     void *ptr_new;
  75
  76     ma(&ptr_new,nbytes_new);
  77
  78     if (nbytes_new > 0 && ptr_new == NULL)
  79     {
  80         gmx_fatal(FARGS, "Allocation of %d bytes failed", nbytes_new);
  81     }
  82
  83     if (nbytes_copy > 0)
  84     {
  85         if (nbytes_new < nbytes_copy)
  86         {
  87             gmx_incons("In nbnxn_realloc_void: new size less than copy size");
  88         }
  89         memcpy(ptr_new,*ptr,nbytes_copy);
  90     }
  91     if (*ptr != NULL)
  92     {
  93         mf(*ptr);
  94     }
  95     *ptr = ptr_new;
  96 }
  97
  98 /* Reallocate the nbnxn_atomdata_t for a size of n atoms */
  99 void nbnxn_atomdata_realloc(nbnxn_atomdata_t *nbat,int n)
 100 {
 101     int t;
 102
 103     nbnxn_realloc_void((void **)&nbat->type,
 104                        nbat->natoms*sizeof(*nbat->type),
 105                        n*sizeof(*nbat->type),
 106                        nbat->alloc,nbat->free);
 107     nbnxn_realloc_void((void **)&nbat->lj_comb,
 108                        nbat->natoms*2*sizeof(*nbat->lj_comb),
 109                        n*2*sizeof(*nbat->lj_comb),
 110                        nbat->alloc,nbat->free);
 111     if (nbat->XFormat != nbatXYZQ)
 112     {
 113         nbnxn_realloc_void((void **)&nbat->q,
 114                            nbat->natoms*sizeof(*nbat->q),
 115                            n*sizeof(*nbat->q),
 116                            nbat->alloc,nbat->free);
 117     }
 118     if (nbat->nenergrp > 1)
 119     {
 120         nbnxn_realloc_void((void **)&nbat->energrp,
 121                            nbat->natoms/nbat->na_c*sizeof(*nbat->energrp),
 122                            n/nbat->na_c*sizeof(*nbat->energrp),
 123                            nbat->alloc,nbat->free);
 124     }
 125     nbnxn_realloc_void((void **)&nbat->x,
 126                        nbat->natoms*nbat->xstride*sizeof(*nbat->x),
 127                        n*nbat->xstride*sizeof(*nbat->x),
 128                        nbat->alloc,nbat->free);
 129     for(t=0; t<nbat->nout; t++)
 130     {
 131         /* Allocate one element extra for possible signaling with CUDA */
 132         nbnxn_realloc_void((void **)&nbat->out[t].f,
 133                            nbat->natoms*nbat->fstride*sizeof(*nbat->out[t].f),
 134                            n*nbat->fstride*sizeof(*nbat->out[t].f),
 135                            nbat->alloc,nbat->free);
 136     }
 137     nbat->nalloc = n;
 138 }
 139
 140 /* Initializes an nbnxn_atomdata_output_t data structure */
 141 static void nbnxn_atomdata_output_init(nbnxn_atomdata_output_t *out,
 142                                        int nb_kernel_type,
 143                                        int nenergrp,int stride,
 144                                        nbnxn_alloc_t *ma)
 145 {
 146     int cj_size;
 147
 148     out->f = NULL;
 149     ma((void **)&out->fshift,SHIFTS*DIM*sizeof(*out->fshift));
 150     out->nV = nenergrp*nenergrp;
 151     ma((void **)&out->Vvdw,out->nV*sizeof(*out->Vvdw));
 152     ma((void **)&out->Vc  ,out->nV*sizeof(*out->Vc  ));
 153
 154     if (nb_kernel_type == nbk4xN_X86_SIMD128 ||
 155         nb_kernel_type == nbk4xN_X86_SIMD256)
 156     {
 157         cj_size = nbnxn_kernel_to_cj_size(nb_kernel_type);
 158         out->nVS = nenergrp*nenergrp*stride*(cj_size>>1)*cj_size;
 159         ma((void **)&out->VSvdw,out->nVS*sizeof(*out->VSvdw));
 160         ma((void **)&out->VSc  ,out->nVS*sizeof(*out->VSc  ));
 161     }
 162     else
 163     {
 164         out->nVS = 0;
 165     }
 166 }
 167
 168 static void copy_int_to_nbat_int(const int *a,int na,int na_round,
 169                                  const int *in,int fill,int *innb)
 170 {
 171     int i,j;
 172
 173     j = 0;
 174     for(i=0; i<na; i++)
 175     {
 176         innb[j++] = in[a[i]];
 177     }
 178     /* Complete the partially filled last cell with fill */
 179     for(; i<na_round; i++)
 180     {
 181         innb[j++] = fill;
 182     }
 183 }
 184
 185 static void clear_nbat_real(int na,int nbatFormat,real *xnb,int a0)
 186 {
 187     int a,d,j,c;
 188
 189     switch (nbatFormat)
 190     {
 191     case nbatXYZ:
 192         for(a=0; a<na; a++)
 193         {
 194             for(d=0; d<DIM; d++)
 195             {
 196                 xnb[(a0+a)*STRIDE_XYZ+d] = 0;
 197             }
 198         }
 199         break;
 200     case nbatXYZQ:
 201         for(a=0; a<na; a++)
 202         {
 203             for(d=0; d<DIM; d++)
 204             {
 205                 xnb[(a0+a)*STRIDE_XYZQ+d] = 0;
 206             }
 207         }
 208         break;
 209     case nbatX4:
 210         j = X4_IND_A(a0);
 211         c = a0 & (PACK_X4-1);
 212         for(a=0; a<na; a++)
 213         {
 214             xnb[j+XX*PACK_X4] = 0;
 215             xnb[j+YY*PACK_X4] = 0;
 216             xnb[j+ZZ*PACK_X4] = 0;
 217             j++;
 218             c++;
 219             if (c == PACK_X4)
 220             {
 221                 j += (DIM-1)*PACK_X4;
 222                 c  = 0;
 223             }
 224         }
 225         break;
 226     case nbatX8:
 227         j = X8_IND_A(a0);
 228         c = a0 & (PACK_X8-1);
 229         for(a=0; a<na; a++)
 230         {
 231             xnb[j+XX*PACK_X8] = 0;
 232             xnb[j+YY*PACK_X8] = 0;
 233             xnb[j+ZZ*PACK_X8] = 0;
 234             j++;
 235             c++;
 236             if (c == PACK_X8)
 237             {
 238                 j += (DIM-1)*PACK_X8;
 239                 c  = 0;
 240             }
 241         }
 242         break;
 243     }
 244 }
 245
 246 void copy_rvec_to_nbat_real(const int *a,int na,int na_round,
 247                             rvec *x,int nbatFormat,real *xnb,int a0,
 248                             int cx,int cy,int cz)
 249 {
 250     int i,j,c;
 251
 252 /* We might need to place filler particles to fill up the cell to na_round.
 253  * The coefficients (LJ and q) for such particles are zero.
 254  * But we might still get NaN as 0*NaN when distances are too small.
 255  * We hope that -107 nm is far away enough from to zero
 256  * to avoid accidental short distances to particles shifted down for pbc.
 257  */
 258 #define NBAT_FAR_AWAY 107
 259
 260     switch (nbatFormat)
 261     {
 262     case nbatXYZ:
 263         j = a0*STRIDE_XYZ;
 264         for(i=0; i<na; i++)
 265         {
 266             xnb[j++] = x[a[i]][XX];
 267             xnb[j++] = x[a[i]][YY];
 268             xnb[j++] = x[a[i]][ZZ];
 269         }
 270         /* Complete the partially filled last cell with copies of the last element.
 271          * This simplifies the bounding box calculation and avoid
 272          * numerical issues with atoms that are coincidentally close.
 273          */
 274         for(; i<na_round; i++)
 275         {
 276             xnb[j++] = -NBAT_FAR_AWAY*(1 + cx);
 277             xnb[j++] = -NBAT_FAR_AWAY*(1 + cy);
 278             xnb[j++] = -NBAT_FAR_AWAY*(1 + cz + i);
 279         }
 280         break;
 281     case nbatXYZQ:
 282         j = a0*STRIDE_XYZQ;
 283         for(i=0; i<na; i++)
 284         {
 285             xnb[j++] = x[a[i]][XX];
 286             xnb[j++] = x[a[i]][YY];
 287             xnb[j++] = x[a[i]][ZZ];
 288             j++;
 289         }
 290         /* Complete the partially filled last cell with particles far apart */
 291         for(; i<na_round; i++)
 292         {
 293             xnb[j++] = -NBAT_FAR_AWAY*(1 + cx);
 294             xnb[j++] = -NBAT_FAR_AWAY*(1 + cy);
 295             xnb[j++] = -NBAT_FAR_AWAY*(1 + cz + i);
 296             j++;
 297         }
 298         break;
 299     case nbatX4:
 300         j = X4_IND_A(a0);
 301         c = a0 & (PACK_X4-1);
 302         for(i=0; i<na; i++)
 303         {
 304             xnb[j+XX*PACK_X4] = x[a[i]][XX];
 305             xnb[j+YY*PACK_X4] = x[a[i]][YY];
 306             xnb[j+ZZ*PACK_X4] = x[a[i]][ZZ];
 307             j++;
 308             c++;
 309             if (c == PACK_X4)
 310             {
 311                 j += (DIM-1)*PACK_X4;
 312                 c  = 0;
 313             }
 314         }
 315         /* Complete the partially filled last cell with particles far apart */
 316         for(; i<na_round; i++)
 317         {
 318             xnb[j+XX*PACK_X4] = -NBAT_FAR_AWAY*(1 + cx);
 319             xnb[j+YY*PACK_X4] = -NBAT_FAR_AWAY*(1 + cy);
 320             xnb[j+ZZ*PACK_X4] = -NBAT_FAR_AWAY*(1 + cz + i);
 321             j++;
 322             c++;
 323             if (c == PACK_X4)
 324             {
 325                 j += (DIM-1)*PACK_X4;
 326                 c  = 0;
 327             }
 328         }
 329         break;
 330     case nbatX8:
 331         j = X8_IND_A(a0);
 332         c = a0 & (PACK_X8 - 1);
 333         for(i=0; i<na; i++)
 334         {
 335             xnb[j+XX*PACK_X8] = x[a[i]][XX];
 336             xnb[j+YY*PACK_X8] = x[a[i]][YY];
 337             xnb[j+ZZ*PACK_X8] = x[a[i]][ZZ];
 338             j++;
 339             c++;
 340             if (c == PACK_X8)
 341             {
 342                 j += (DIM-1)*PACK_X8;
 343                 c  = 0;
 344             }
 345         }
 346         /* Complete the partially filled last cell with particles far apart */
 347         for(; i<na_round; i++)
 348         {
 349             xnb[j+XX*PACK_X8] = -NBAT_FAR_AWAY*(1 + cx);
 350             xnb[j+YY*PACK_X8] = -NBAT_FAR_AWAY*(1 + cy);
 351             xnb[j+ZZ*PACK_X8] = -NBAT_FAR_AWAY*(1 + cz + i);
 352             j++;
 353             c++;
 354             if (c == PACK_X8)
 355             {
 356                 j += (DIM-1)*PACK_X8;
 357                 c  = 0;
 358             }
 359         }
 360         break;
 361     default:
 362         gmx_incons("Unsupported stride");
 363     }
 364 }
 365
 366 /* Determines the combination rule (or none) to be used, stores it,
 367  * and sets the LJ parameters required with the rule.
 368  */
 369 static void set_combination_rule_data(nbnxn_atomdata_t *nbat)
 370 {
 371     int  nt,i,j;
 372     real c6,c12;
 373
 374     nt = nbat->ntype;
 375
 376     switch (nbat->comb_rule)
 377     {
 378     case  ljcrGEOM:
 379         nbat->comb_rule = ljcrGEOM;
 380
 381         for(i=0; i<nt; i++)
 382         {
 383             /* Copy the diagonal from the nbfp matrix */
 384             nbat->nbfp_comb[i*2  ] = sqrt(nbat->nbfp[(i*nt+i)*2  ]);
 385             nbat->nbfp_comb[i*2+1] = sqrt(nbat->nbfp[(i*nt+i)*2+1]);
 386         }
 387         break;
 388     case ljcrLB:
 389         for(i=0; i<nt; i++)
 390         {
 391             /* Get 6*C6 and 12*C12 from the diagonal of the nbfp matrix */
 392             c6  = nbat->nbfp[(i*nt+i)*2  ];
 393             c12 = nbat->nbfp[(i*nt+i)*2+1];
 394             if (c6 > 0 && c12 > 0)
 395             {
 396                 /* We store 0.5*2^1/6*sigma and sqrt(4*3*eps),
 397                  * so we get 6*C6 and 12*C12 after combining.
 398                  */
 399                 nbat->nbfp_comb[i*2  ] = 0.5*pow(c12/c6,1.0/6.0);
 400                 nbat->nbfp_comb[i*2+1] = sqrt(c6*c6/c12);
 401             }
 402             else
 403             {
 404                 nbat->nbfp_comb[i*2  ] = 0;
 405                 nbat->nbfp_comb[i*2+1] = 0;
 406             }
 407         }
 408         break;
 409     case ljcrNONE:
 410         /* In nbfp_s4 we use a stride of 4 for storing two parameters */
 411         nbat->alloc((void **)&nbat->nbfp_s4,nt*nt*4*sizeof(*nbat->nbfp_s4));
 412         for(i=0; i<nt; i++)
 413         {
 414             for(j=0; j<nt; j++)
 415             {
 416                 nbat->nbfp_s4[(i*nt+j)*4+0] = nbat->nbfp[(i*nt+j)*2+0];
 417                 nbat->nbfp_s4[(i*nt+j)*4+1] = nbat->nbfp[(i*nt+j)*2+1];
 418                 nbat->nbfp_s4[(i*nt+j)*4+2] = 0;
 419                 nbat->nbfp_s4[(i*nt+j)*4+3] = 0;
 420             }
 421         }
 422         break;
 423     default:
 424         gmx_incons("Unknown combination rule");
 425         break;
 426     }
 427 }
 428
 429 /* Initializes an nbnxn_atomdata_t data structure */
 430 void nbnxn_atomdata_init(FILE *fp,
 431                          nbnxn_atomdata_t *nbat,
 432                          int nb_kernel_type,
 433                          int ntype,const real *nbfp,
 434                          int n_energygroups,
 435                          int nout,
 436                          nbnxn_alloc_t *alloc,
 437                          nbnxn_free_t  *free)
 438 {
 439     int  i,j;
 440     real c6,c12,tol;
 441     char *ptr;
 442     gmx_bool simple,bCombGeom,bCombLB;
 443
 444     if (alloc == NULL)
 445     {
 446         nbat->alloc = nbnxn_alloc_aligned;
 447     }
 448     else
 449     {
 450         nbat->alloc = alloc;
 451     }
 452     if (free == NULL)
 453     {
 454         nbat->free = nbnxn_free_aligned;
 455     }
 456     else
 457     {
 458         nbat->free = free;
 459     }
 460
 461     if (debug)
 462     {
 463         fprintf(debug,"There are %d atom types in the system, adding one for nbnxn_atomdata_t\n",ntype);
 464     }
 465     nbat->ntype = ntype + 1;
 466     nbat->alloc((void **)&nbat->nbfp,
 467                 nbat->ntype*nbat->ntype*2*sizeof(*nbat->nbfp));
 468     nbat->alloc((void **)&nbat->nbfp_comb,nbat->ntype*2*sizeof(*nbat->nbfp_comb));
 469
 470     /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
 471      * force-field floating point parameters.
 472      */
 473     tol = 1e-5;
 474     ptr = getenv("GMX_LJCOMB_TOL");
 475     if (ptr != NULL)
 476     {
 477         double dbl;
 478
 479         sscanf(ptr,"%lf",&dbl);
 480         tol = dbl;
 481     }
 482     bCombGeom = TRUE;
 483     bCombLB   = TRUE;
 484
 485     /* Temporarily fill nbat->nbfp_comb with sigma and epsilon
 486      * to check for the LB rule.
 487      */
 488     for(i=0; i<ntype; i++)
 489     {
 490         c6  = nbfp[(i*ntype+i)*2  ]/6.0;
 491         c12 = nbfp[(i*ntype+i)*2+1]/12.0;
 492         if (c6 > 0 && c12 > 0)
 493         {
 494             nbat->nbfp_comb[i*2  ] = pow(c12/c6,1.0/6.0);
 495             nbat->nbfp_comb[i*2+1] = 0.25*c6*c6/c12;
 496         }
 497         else if (c6 == 0 && c12 == 0)
 498         {
 499             nbat->nbfp_comb[i*2  ] = 0;
 500             nbat->nbfp_comb[i*2+1] = 0;
 501         }
 502         else
 503         {
 504             /* Can not use LB rule with only dispersion or repulsion */
 505             bCombLB = FALSE;
 506         }
 507     }
 508
 509     for(i=0; i<nbat->ntype; i++)
 510     {
 511         for(j=0; j<nbat->ntype; j++)
 512         {
 513             if (i < ntype && j < ntype)
 514             {
 515                 /* fr->nbfp has been updated, so that array too now stores c6/c12 including
 516                  * the 6.0/12.0 prefactors to save 2 flops in the most common case (force-only).
 517                  */
 518                 c6  = nbfp[(i*ntype+j)*2  ];
 519                 c12 = nbfp[(i*ntype+j)*2+1];
 520                 nbat->nbfp[(i*nbat->ntype+j)*2  ] = c6;
 521                 nbat->nbfp[(i*nbat->ntype+j)*2+1] = c12;
 522
 523                 /* Compare 6*C6 and 12*C12 for geometric cobination rule */
 524                 bCombGeom = bCombGeom &&
 525                     gmx_within_tol(c6*c6  ,nbfp[(i*ntype+i)*2  ]*nbfp[(j*ntype+j)*2  ],tol) &&
 526                     gmx_within_tol(c12*c12,nbfp[(i*ntype+i)*2+1]*nbfp[(j*ntype+j)*2+1],tol);
 527
 528                 /* Compare C6 and C12 for Lorentz-Berthelot combination rule */
 529                 c6  /= 6.0;
 530                 c12 /= 12.0;
 531                 bCombLB = bCombLB &&
 532                     ((c6 == 0 && c12 == 0 &&
 533                       (nbat->nbfp_comb[i*2+1] == 0 || nbat->nbfp_comb[j*2+1] == 0)) ||
 534                      (c6 > 0 && c12 > 0 &&
 535                       gmx_within_tol(pow(c12/c6,1.0/6.0),0.5*(nbat->nbfp_comb[i*2]+nbat->nbfp_comb[j*2]),tol) &&
 536                       gmx_within_tol(0.25*c6*c6/c12,sqrt(nbat->nbfp_comb[i*2+1]*nbat->nbfp_comb[j*2+1]),tol)));
 537             }
 538             else
 539             {
 540                 /* Add zero parameters for the additional dummy atom type */
 541                 nbat->nbfp[(i*nbat->ntype+j)*2  ] = 0;
 542                 nbat->nbfp[(i*nbat->ntype+j)*2+1] = 0;
 543             }
 544         }
 545     }
 546     if (debug)
 547     {
 548         fprintf(debug,"Combination rules: geometric %d Lorentz-Berthelot %d\n",
 549                 bCombGeom,bCombLB);
 550     }
 551
 552     simple = nbnxn_kernel_pairlist_simple(nb_kernel_type);
 553
 554     if (simple)
 555     {
 556         /* We prefer the geometic combination rule,
 557          * as that gives a slightly faster kernel than the LB rule.
 558          */
 559         if (bCombGeom)
 560         {
 561             nbat->comb_rule = ljcrGEOM;
 562         }
 563         else if (bCombLB)
 564         {
 565             nbat->comb_rule = ljcrLB;
 566         }
 567         else
 568         {
 569             nbat->comb_rule = ljcrNONE;
 570
 571             nbat->free(nbat->nbfp_comb);
 572         }
 573
 574         if (fp)
 575         {
 576             if (nbat->comb_rule == ljcrNONE)
 577             {
 578                 fprintf(fp,"Using full Lennard-Jones parameter combination matrix\n\n");
 579             }
 580             else
 581             {
 582                 fprintf(fp,"Using %s Lennard-Jones combination rule\n\n",
 583                         nbat->comb_rule==ljcrGEOM ? "geometric" : "Lorentz-Berthelot");
 584             }
 585         }
 586
 587         set_combination_rule_data(nbat);
 588     }
 589     else
 590     {
 591         nbat->comb_rule = ljcrNONE;
 592
 593         nbat->free(nbat->nbfp_comb);
 594     }
 595
 596     nbat->natoms  = 0;
 597     nbat->type    = NULL;
 598     nbat->lj_comb = NULL;
 599     if (simple)
 600     {
 601         switch (nb_kernel_type)
 602         {
 603         case nbk4xN_X86_SIMD128:
 604             nbat->XFormat = nbatX4;
 605             break;
 606         case nbk4xN_X86_SIMD256:
 607 #ifndef GMX_DOUBLE
 608             nbat->XFormat = nbatX8;
 609 #else
 610             nbat->XFormat = nbatX4;
 611 #endif
 612             break;
 613         default:
 614             nbat->XFormat = nbatXYZ;
 615             break;
 616         }
 617
 618         nbat->FFormat = nbat->XFormat;
 619     }
 620     else
 621     {
 622         nbat->XFormat = nbatXYZQ;
 623         nbat->FFormat = nbatXYZ;
 624     }
 625     nbat->q       = NULL;
 626     nbat->nenergrp = n_energygroups;
 627     if (!simple)
 628     {
 629         /* Energy groups not supported yet for super-sub lists */
 630         if (n_energygroups > 1 && fp != NULL)
 631         {
 632             fprintf(fp,"\nNOTE: With GPUs, reporting energy group contributions is not supported\n\n");
 633         }
 634         nbat->nenergrp = 1;
 635     }
 636     /* Temporary storage goes as #grp^3*simd_width^2/2, so limit to 64 */
 637     if (nbat->nenergrp > 64)
 638     {
 639         gmx_fatal(FARGS,"With NxN kernels not more than 64 energy groups are supported\n");
 640     }
 641     nbat->neg_2log = 1;
 642     while (nbat->nenergrp > (1<<nbat->neg_2log))
 643     {
 644         nbat->neg_2log++;
 645     }
 646     nbat->energrp = NULL;
 647     nbat->alloc((void **)&nbat->shift_vec,SHIFTS*sizeof(*nbat->shift_vec));
 648     nbat->xstride = (nbat->XFormat == nbatXYZQ ? STRIDE_XYZQ : DIM);
 649     nbat->fstride = (nbat->FFormat == nbatXYZQ ? STRIDE_XYZQ : DIM);
 650     nbat->x       = NULL;
 651     nbat->nout    = nout;
 652     snew(nbat->out,nbat->nout);
 653     nbat->nalloc  = 0;
 654     for(i=0; i<nbat->nout; i++)
 655     {
 656         nbnxn_atomdata_output_init(&nbat->out[i],
 657                                    nb_kernel_type,
 658                                    nbat->nenergrp,1<<nbat->neg_2log,
 659                                    nbat->alloc);
 660     }
 661 }
 662
 663 static void copy_lj_to_nbat_lj_comb_x4(const real *ljparam_type,
 664                                        const int *type,int na,
 665                                        real *ljparam_at)
 666 {
 667     int is,k,i;
 668
 669     /* The LJ params follow the combination rule:
 670      * copy the params for the type array to the atom array.
 671      */
 672     for(is=0; is<na; is+=PACK_X4)
 673     {
 674         for(k=0; k<PACK_X4; k++)
 675         {
 676             i = is + k;
 677             ljparam_at[is*2        +k] = ljparam_type[type[i]*2  ];
 678             ljparam_at[is*2+PACK_X4+k] = ljparam_type[type[i]*2+1];
 679         }
 680     }
 681 }
 682
 683 static void copy_lj_to_nbat_lj_comb_x8(const real *ljparam_type,
 684                                        const int *type,int na,
 685                                        real *ljparam_at)
 686 {
 687     int is,k,i;
 688
 689     /* The LJ params follow the combination rule:
 690      * copy the params for the type array to the atom array.
 691      */
 692     for(is=0; is<na; is+=PACK_X8)
 693     {
 694         for(k=0; k<PACK_X8; k++)
 695         {
 696             i = is + k;
 697             ljparam_at[is*2        +k] = ljparam_type[type[i]*2  ];
 698             ljparam_at[is*2+PACK_X8+k] = ljparam_type[type[i]*2+1];
 699         }
 700     }
 701 }
 702
 703 /* Sets the atom type and LJ data in nbnxn_atomdata_t */
 704 static void nbnxn_atomdata_set_atomtypes(nbnxn_atomdata_t *nbat,
 705                                          int ngrid,
 706                                          const nbnxn_search_t nbs,
 707                                          const int *type)
 708 {
 709     int g,i,ncz,ash;
 710     const nbnxn_grid_t *grid;
 711
 712     for(g=0; g<ngrid; g++)
 713     {
 714         grid = &nbs->grid[g];
 715
 716         /* Loop over all columns and copy and fill */
 717         for(i=0; i<grid->ncx*grid->ncy; i++)
 718         {
 719             ncz = grid->cxy_ind[i+1] - grid->cxy_ind[i];
 720             ash = (grid->cell0 + grid->cxy_ind[i])*grid->na_sc;
 721
 722             copy_int_to_nbat_int(nbs->a+ash,grid->cxy_na[i],ncz*grid->na_sc,
 723                                  type,nbat->ntype-1,nbat->type+ash);
 724
 725             if (nbat->comb_rule != ljcrNONE)
 726             {
 727                 if (nbat->XFormat == nbatX4)
 728                 {
 729                     copy_lj_to_nbat_lj_comb_x4(nbat->nbfp_comb,
 730                                                nbat->type+ash,ncz*grid->na_sc,
 731                                                nbat->lj_comb+ash*2);
 732                 }
 733                 else if (nbat->XFormat == nbatX8)
 734                 {
 735                     copy_lj_to_nbat_lj_comb_x8(nbat->nbfp_comb,
 736                                                nbat->type+ash,ncz*grid->na_sc,
 737                                                nbat->lj_comb+ash*2);
 738                 }
 739             }
 740         }
 741     }
 742 }
 743
 744 /* Sets the charges in nbnxn_atomdata_t *nbat */
 745 static void nbnxn_atomdata_set_charges(nbnxn_atomdata_t *nbat,
 746                                        int ngrid,
 747                                        const nbnxn_search_t nbs,
 748                                        const real *charge)
 749 {
 750     int  g,cxy,ncz,ash,na,na_round,i,j;
 751     real *q;
 752     const nbnxn_grid_t *grid;
 753
 754     for(g=0; g<ngrid; g++)
 755     {
 756         grid = &nbs->grid[g];
 757
 758         /* Loop over all columns and copy and fill */
 759         for(cxy=0; cxy<grid->ncx*grid->ncy; cxy++)
 760         {
 761             ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
 762             na  = grid->cxy_na[cxy];
 763             na_round = (grid->cxy_ind[cxy+1] - grid->cxy_ind[cxy])*grid->na_sc;
 764
 765             if (nbat->XFormat == nbatXYZQ)
 766             {
 767                 q = nbat->x + ash*STRIDE_XYZQ + ZZ + 1;
 768                 for(i=0; i<na; i++)
 769                 {
 770                     *q = charge[nbs->a[ash+i]];
 771                     q += STRIDE_XYZQ;
 772                 }
 773                 /* Complete the partially filled last cell with zeros */
 774                 for(; i<na_round; i++)
 775                 {
 776                     *q = 0;
 777                     q += STRIDE_XYZQ;
 778                 }
 779             }
 780             else
 781             {
 782                 q = nbat->q + ash;
 783                 for(i=0; i<na; i++)
 784                 {
 785                     *q = charge[nbs->a[ash+i]];
 786                     q++;
 787                 }
 788                 /* Complete the partially filled last cell with zeros */
 789                 for(; i<na_round; i++)
 790                 {
 791                     *q = 0;
 792                     q++;
 793                 }
 794             }
 795         }
 796     }
 797 }
 798
 799 /* Copies the energy group indices to a reordered and packed array */
 800 static void copy_egp_to_nbat_egps(const int *a,int na,int na_round,
 801                                   int na_c,int bit_shift,
 802                                   const int *in,int *innb)
 803 {
 804     int i,j,sa,at;
 805     int comb;
 806
 807     j = 0;
 808     for(i=0; i<na; i+=na_c)
 809     {
 810         /* Store na_c energy group numbers into one int */
 811         comb = 0;
 812         for(sa=0; sa<na_c; sa++)
 813         {
 814             at = a[i+sa];
 815             if (at >= 0)
 816             {
 817                 comb |= (GET_CGINFO_GID(in[at]) << (sa*bit_shift));
 818             }
 819         }
 820         innb[j++] = comb;
 821     }
 822     /* Complete the partially filled last cell with fill */
 823     for(; i<na_round; i+=na_c)
 824     {
 825         innb[j++] = 0;
 826     }
 827 }
 828
 829 /* Set the energy group indices for atoms in nbnxn_atomdata_t */
 830 static void nbnxn_atomdata_set_energygroups(nbnxn_atomdata_t *nbat,
 831                                             int ngrid,
 832                                             const nbnxn_search_t nbs,
 833                                             const int *atinfo)
 834 {
 835     int g,i,ncz,ash;
 836     const nbnxn_grid_t *grid;
 837
 838     for(g=0; g<ngrid; g++)
 839     {
 840         grid = &nbs->grid[g];
 841
 842         /* Loop over all columns and copy and fill */
 843         for(i=0; i<grid->ncx*grid->ncy; i++)
 844         {
 845             ncz = grid->cxy_ind[i+1] - grid->cxy_ind[i];
 846             ash = (grid->cell0 + grid->cxy_ind[i])*grid->na_sc;
 847
 848             copy_egp_to_nbat_egps(nbs->a+ash,grid->cxy_na[i],ncz*grid->na_sc,
 849                                   nbat->na_c,nbat->neg_2log,
 850                                   atinfo,nbat->energrp+(ash>>grid->na_c_2log));
 851         }
 852     }
 853 }
 854
 855 /* Sets all required atom parameter data in nbnxn_atomdata_t */
 856 void nbnxn_atomdata_set(nbnxn_atomdata_t *nbat,
 857                         int locality,
 858                         const nbnxn_search_t nbs,
 859                         const t_mdatoms *mdatoms,
 860                         const int *atinfo)
 861 {
 862     int ngrid;
 863
 864     if (locality == eatLocal)
 865     {
 866         ngrid = 1;
 867     }
 868     else
 869     {
 870         ngrid = nbs->ngrid;
 871     }
 872
 873     nbnxn_atomdata_set_atomtypes(nbat,ngrid,nbs,mdatoms->typeA);
 874
 875     nbnxn_atomdata_set_charges(nbat,ngrid,nbs,mdatoms->chargeA);
 876
 877     if (nbat->nenergrp > 1)
 878     {
 879         nbnxn_atomdata_set_energygroups(nbat,ngrid,nbs,atinfo);
 880     }
 881 }
 882
 883 /* Copies the shift vector array to nbnxn_atomdata_t */
 884 void nbnxn_atomdata_copy_shiftvec(gmx_bool bDynamicBox,
 885                                    rvec *shift_vec,
 886                                    nbnxn_atomdata_t *nbat)
 887 {
 888     int i;
 889
 890     nbat->bDynamicBox = bDynamicBox;
 891     for(i=0; i<SHIFTS; i++)
 892     {
 893         copy_rvec(shift_vec[i],nbat->shift_vec[i]);
 894     }
 895 }
 896
 897 /* Copies (and reorders) the coordinates to nbnxn_atomdata_t */
 898 void nbnxn_atomdata_copy_x_to_nbat_x(const nbnxn_search_t nbs,
 899                                       int locality,
 900                                       gmx_bool FillLocal,
 901                                       rvec *x,
 902                                       nbnxn_atomdata_t *nbat)
 903 {
 904     int g0=0,g1=0;
 905     int nth,th;
 906
 907     switch (locality)
 908     {
 909     case eatAll:
 910         g0 = 0;
 911         g1 = nbs->ngrid;
 912         break;
 913     case eatLocal:
 914         g0 = 0;
 915         g1 = 1;
 916         break;
 917     case eatNonlocal:
 918         g0 = 1;
 919         g1 = nbs->ngrid;
 920         break;
 921     }
 922
 923     if (FillLocal)
 924     {
 925         nbat->natoms_local = nbs->grid[0].nc*nbs->grid[0].na_sc;
 926     }
 927
 928     nth = gmx_omp_nthreads_get(emntPairsearch);
 929
 930 #pragma omp parallel for num_threads(nth) schedule(static)
 931     for(th=0; th<nth; th++)
 932     {
 933         int g;
 934
 935         for(g=g0; g<g1; g++)
 936         {
 937             const nbnxn_grid_t *grid;
 938             int cxy0,cxy1,cxy;
 939
 940             grid = &nbs->grid[g];
 941
 942             cxy0 = (grid->ncx*grid->ncy* th   +nth-1)/nth;
 943             cxy1 = (grid->ncx*grid->ncy*(th+1)+nth-1)/nth;
 944
 945             for(cxy=cxy0; cxy<cxy1; cxy++)
 946             {
 947                 int na,ash,na_fill;
 948
 949                 na  = grid->cxy_na[cxy];
 950                 ash = (grid->cell0 + grid->cxy_ind[cxy])*grid->na_sc;
 951
 952                 if (g == 0 && FillLocal)
 953                 {
 954                     na_fill =
 955                         (grid->cxy_ind[cxy+1] - grid->cxy_ind[cxy])*grid->na_sc;
 956                 }
 957                 else
 958                 {
 959                     /* We fill only the real particle locations.
 960                      * We assume the filling entries at the end have been
 961                      * properly set before during ns.
 962                      */
 963                     na_fill = na;
 964                 }
 965                 copy_rvec_to_nbat_real(nbs->a+ash,na,na_fill,x,
 966                                        nbat->XFormat,nbat->x,ash,
 967                                        0,0,0);
 968             }
 969         }
 970     }
 971 }
 972
 973 static void
 974 nbnxn_atomdata_reduce_reals(real * gmx_restrict dest,
 975                             real ** gmx_restrict src,
 976                             int nsrc,
 977                             int i0, int i1)
 978 {
 979     int i,s;
 980
 981     for(i=i0; i<i1; i++)
 982     {
 983         for(s=0; s<nsrc; s++)
 984         {
 985             dest[i] += src[s][i];
 986         }
 987     }
 988 }
 989
 990 static void
 991 nbnxn_atomdata_reduce_reals_x86_simd(real * gmx_restrict dest,
 992                                      real ** gmx_restrict src,
 993                                      int nsrc,
 994                                      int i0, int i1)
 995 {
 996 #ifdef NBNXN_SEARCH_SSE
 997 /* We can use AVX256 here, but not when AVX128 kernels are selected.
 998  * As this reduction is not faster with AVX256 anyway, we use 128-bit SIMD.
 999  */
1000 #define GMX_MM128_HERE
1001 #include "gmx_x86_simd_macros.h"
1002
1003     int       i,s;
1004     gmx_mm_pr dest_SSE,src_SSE;
1005
1006     if ((i0 & (GMX_X86_SIMD_WIDTH_HERE-1)) ||
1007         (i1 & (GMX_X86_SIMD_WIDTH_HERE-1)))
1008     {
1009         gmx_incons("bounds not a multiple of GMX_X86_SIMD_WIDTH_HERE in nbnxn_atomdata_reduce_reals_x86_simd");
1010     }
1011
1012     for(i=i0; i<i1; i+=GMX_X86_SIMD_WIDTH_HERE)
1013     {
1014         dest_SSE = gmx_load_pr(dest+i);
1015         for(s=0; s<nsrc; s++)
1016         {
1017             src_SSE  = gmx_load_pr(src[s]+i);
1018             dest_SSE = gmx_add_pr(dest_SSE,src_SSE);
1019         }
1020         gmx_store_pr(dest+i,dest_SSE);
1021     }
1022
1023 #undef GMX_MM128_HERE
1024 #undef GMX_MM256_HERE
1025 #endif
1026 }
1027
1028 /* Add part of the force array(s) from nbnxn_atomdata_t to f */
1029 static void
1030 nbnxn_atomdata_add_nbat_f_to_f_part(const nbnxn_search_t nbs,
1031                                     const nbnxn_atomdata_t *nbat,
1032                                     nbnxn_atomdata_output_t *out,
1033                                     int nfa,
1034                                     int a0,int a1,
1035                                     rvec *f)
1036 {
1037     int  a,i,fa;
1038     const int  *cell;
1039     const real *fnb;
1040
1041     cell = nbs->cell;
1042
1043     /* Loop over all columns and copy and fill */
1044     switch (nbat->FFormat)
1045     {
1046     case nbatXYZ:
1047     case nbatXYZQ:
1048         if (nfa == 1)
1049         {
1050             fnb = out[0].f;
1051
1052             for(a=a0; a<a1; a++)
1053             {
1054                 i = cell[a]*nbat->fstride;
1055
1056                 f[a][XX] += fnb[i];
1057                 f[a][YY] += fnb[i+1];
1058                 f[a][ZZ] += fnb[i+2];
1059             }
1060         }
1061         else
1062         {
1063             for(a=a0; a<a1; a++)
1064             {
1065                 i = cell[a]*nbat->fstride;
1066
1067                 for(fa=0; fa<nfa; fa++)
1068                 {
1069                     f[a][XX] += out[fa].f[i];
1070                     f[a][YY] += out[fa].f[i+1];
1071                     f[a][ZZ] += out[fa].f[i+2];
1072                 }
1073             }
1074         }
1075         break;
1076     case nbatX4:
1077         if (nfa == 1)
1078         {
1079             fnb = out[0].f;
1080
1081             for(a=a0; a<a1; a++)
1082             {
1083                 i = X4_IND_A(cell[a]);
1084
1085                 f[a][XX] += fnb[i+XX*PACK_X4];
1086                 f[a][YY] += fnb[i+YY*PACK_X4];
1087                 f[a][ZZ] += fnb[i+ZZ*PACK_X4];
1088             }
1089         }
1090         else
1091         {
1092             for(a=a0; a<a1; a++)
1093             {
1094                 i = X4_IND_A(cell[a]);
1095
1096                 for(fa=0; fa<nfa; fa++)
1097                 {
1098                     f[a][XX] += out[fa].f[i+XX*PACK_X4];
1099                     f[a][YY] += out[fa].f[i+YY*PACK_X4];
1100                     f[a][ZZ] += out[fa].f[i+ZZ*PACK_X4];
1101                 }
1102             }
1103         }
1104         break;
1105     case nbatX8:
1106         if (nfa == 1)
1107         {
1108             fnb = out[0].f;
1109
1110             for(a=a0; a<a1; a++)
1111             {
1112                 i = X8_IND_A(cell[a]);
1113
1114                 f[a][XX] += fnb[i+XX*PACK_X8];
1115                 f[a][YY] += fnb[i+YY*PACK_X8];
1116                 f[a][ZZ] += fnb[i+ZZ*PACK_X8];
1117             }
1118         }
1119         else
1120         {
1121             for(a=a0; a<a1; a++)
1122             {
1123                 i = X8_IND_A(cell[a]);
1124
1125                 for(fa=0; fa<nfa; fa++)
1126                 {
1127                     f[a][XX] += out[fa].f[i+XX*PACK_X8];
1128                     f[a][YY] += out[fa].f[i+YY*PACK_X8];
1129                     f[a][ZZ] += out[fa].f[i+ZZ*PACK_X8];
1130                 }
1131             }
1132         }
1133         break;
1134     }
1135 }
1136
1137 /* Add the force array(s) from nbnxn_atomdata_t to f */
1138 void nbnxn_atomdata_add_nbat_f_to_f(const nbnxn_search_t nbs,
1139                                     int locality,
1140                                     const nbnxn_atomdata_t *nbat,
1141                                     rvec *f)
1142 {
1143     int a0=0,na=0;
1144     int nth,th;
1145     gmx_bool bStreamingReduce;
1146
1147     nbs_cycle_start(&nbs->cc[enbsCCreducef]);
1148
1149     switch (locality)
1150     {
1151     case eatAll:
1152         a0 = 0;
1153         na = nbs->natoms_nonlocal;
1154         break;
1155     case eatLocal:
1156         a0 = 0;
1157         na = nbs->natoms_local;
1158         break;
1159     case eatNonlocal:
1160         a0 = nbs->natoms_local;
1161         na = nbs->natoms_nonlocal - nbs->natoms_local;
1162         break;
1163     }
1164
1165     nth = gmx_omp_nthreads_get(emntNonbonded);
1166
1167     /* Using the two-step streaming reduction is probably always faster */
1168     bStreamingReduce = (nbat->nout > 1);
1169
1170     if (bStreamingReduce)
1171     {
1172         /* Reduce the force thread output buffers into buffer 0, before adding
1173          * them to the, differently ordered, "real" force buffer.
1174          */
1175 #pragma omp parallel for num_threads(nth) schedule(static)
1176         for(th=0; th<nth; th++)
1177         {
1178             int g0,g1,g;
1179
1180             /* For which grids should we reduce the force output? */
1181             g0 = ((locality==eatLocal || locality==eatAll) ? 0 : 1);
1182             g1 = (locality==eatLocal ? 1 : nbs->ngrid);
1183
1184             for(g=g0; g<g1; g++)
1185             {
1186                 nbnxn_grid_t *grid;
1187                 int b0,b1,b;
1188                 int c0,c1,i0,i1;
1189                 int nfptr;
1190                 real *fptr[NBNXN_CELLBLOCK_MAX_THREADS];
1191                 int out;
1192
1193                 grid = &nbs->grid[g];
1194
1195                 /* Calculate the cell-block range for our thread */
1196                 b0 = (grid->cellblock_flags.ncb* th   )/nth;
1197                 b1 = (grid->cellblock_flags.ncb*(th+1))/nth;
1198
1199                 if (grid->cellblock_flags.bUse)
1200                 {
1201                     for(b=b0; b<b1; b++)
1202                     {
1203                         c0 = b*NBNXN_CELLBLOCK_SIZE;
1204                         c1 = min(c0 + NBNXN_CELLBLOCK_SIZE,grid->nc);
1205                         i0 = (grid->cell0 + c0)*grid->na_c*nbat->fstride;
1206                         i1 = (grid->cell0 + c1)*grid->na_c*nbat->fstride;
1207
1208                         nfptr = 0;
1209                         for(out=1; out<nbat->nout; out++)
1210                         {
1211                             if (grid->cellblock_flags.flag[b] & (1U<<out))
1212                             {
1213                                 fptr[nfptr++] = nbat->out[out].f;
1214                             }
1215                         }
1216                         if (nfptr > 0)
1217                         {
1218 #ifdef NBNXN_SEARCH_SSE
1219                             nbnxn_atomdata_reduce_reals_x86_simd
1220 #else
1221                             nbnxn_atomdata_reduce_reals
1222 #endif
1223                                                        (nbat->out[0].f,
1224                                                         fptr,nfptr,
1225                                                         i0,i1);
1226                         }
1227                     }
1228                 }
1229                 else
1230                 {
1231                     c0 = b0*NBNXN_CELLBLOCK_SIZE;
1232                     c1 = min(b1*NBNXN_CELLBLOCK_SIZE,grid->nc);
1233                     i0 = (grid->cell0 + c0)*grid->na_c*nbat->fstride;
1234                     i1 = (grid->cell0 + c1)*grid->na_c*nbat->fstride;
1235
1236                     nfptr = 0;
1237                     for(out=1; out<nbat->nout; out++)
1238                     {
1239                         fptr[nfptr++] = nbat->out[out].f;
1240                     }
1241
1242 #ifdef NBNXN_SEARCH_SSE
1243                     nbnxn_atomdata_reduce_reals_x86_simd
1244 #else
1245                     nbnxn_atomdata_reduce_reals
1246 #endif
1247                                                (nbat->out[0].f,
1248                                                 fptr,nfptr,
1249                                                 i0,i1);
1250                 }
1251             }
1252         }
1253     }
1254
1255 #pragma omp parallel for num_threads(nth) schedule(static)
1256     for(th=0; th<nth; th++)
1257     {
1258         nbnxn_atomdata_add_nbat_f_to_f_part(nbs,nbat,
1259                                             nbat->out,
1260                                             bStreamingReduce ? 1 : nbat->nout,
1261                                             a0+((th+0)*na)/nth,
1262                                             a0+((th+1)*na)/nth,
1263                                             f);
1264     }
1265
1266     nbs_cycle_stop(&nbs->cc[enbsCCreducef]);
1267 }
1268
1269 /* Adds the shift forces from nbnxn_atomdata_t to fshift */
1270 void nbnxn_atomdata_add_nbat_fshift_to_fshift(const nbnxn_atomdata_t *nbat,
1271                                               rvec *fshift)
1272 {
1273     const nbnxn_atomdata_output_t *out;
1274     int  th;
1275     int  s;
1276     rvec sum;
1277
1278     out = nbat->out;
1279
1280     for(s=0; s<SHIFTS; s++)
1281     {
1282         clear_rvec(sum);
1283         for(th=0; th<nbat->nout; th++)
1284         {
1285             sum[XX] += out[th].fshift[s*DIM+XX];
1286             sum[YY] += out[th].fshift[s*DIM+YY];
1287             sum[ZZ] += out[th].fshift[s*DIM+ZZ];
1288         }
1289         rvec_inc(fshift[s],sum);
1290     }
1291 }