src/mdlib/ns.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #ifdef GMX_THREAD_SHM_FDECOMP
  41 #include <pthread.h>
  42 #endif
  43
  44 #include <math.h>
  45 #include <string.h>
  46 #include "sysstuff.h"
  47 #include "smalloc.h"
  48 #include "macros.h"
  49 #include "maths.h"
  50 #include "vec.h"
  51 #include "network.h"
  52 #include "nsgrid.h"
  53 #include "force.h"
  54 #include "nonbonded.h"
  55 #include "ns.h"
  56 #include "pbc.h"
  57 #include "names.h"
  58 #include "gmx_fatal.h"
  59 #include "nrnb.h"
  60 #include "txtdump.h"
  61 #include "mtop_util.h"
  62
  63 #include "domdec.h"
  64
  65
  66 /*
  67  *    E X C L U S I O N   H A N D L I N G
  68  */
  69
  70 #ifdef DEBUG
  71 static void SETEXCL_(t_excl e[],atom_id i,atom_id j)
  72 {   e[j] = e[j] | (1<<i); }
  73 static void RMEXCL_(t_excl e[],atom_id i,atom_id j)
  74 { e[j]=e[j] & ~(1<<i); }
  75 static bool ISEXCL_(t_excl e[],atom_id i,atom_id j)
  76 { return (bool)(e[j] & (1<<i)); }
  77 static bool NOTEXCL_(t_excl e[],atom_id i,atom_id j)
  78 {  return !(ISEXCL(e,i,j)); }
  79 #else
  80 #define SETEXCL(e,i,j) (e)[((atom_id) (j))] |= (1<<((atom_id) (i)))
  81 #define RMEXCL(e,i,j)  (e)[((atom_id) (j))] &= (~(1<<((atom_id) (i))))
  82 #define ISEXCL(e,i,j)  (bool) ((e)[((atom_id) (j))] & (1<<((atom_id) (i))))
  83 #define NOTEXCL(e,i,j) !(ISEXCL(e,i,j))
  84 #endif
  85
  86 /************************************************
  87  *
  88  *  U T I L I T I E S    F O R    N S
  89  *
  90  ************************************************/
  91
  92 static void reallocate_nblist(t_nblist *nl)
  93 {
  94     if (gmx_debug_at)
  95     {
  96         fprintf(debug,"reallocating neigborlist il_code=%d, maxnri=%d\n",
  97                 nl->il_code,nl->maxnri);
  98     }
  99     srenew(nl->iinr,   nl->maxnri);
 100     if (nl->enlist == enlistCG_CG)
 101     {
 102         srenew(nl->iinr_end,nl->maxnri);
 103     }
 104     srenew(nl->gid,    nl->maxnri);
 105     srenew(nl->shift,  nl->maxnri);
 106     srenew(nl->jindex, nl->maxnri+1);
 107 }
 108
 109 /* ivdw/icoul are used to determine the type of interaction, so we
 110  * can set an innerloop index here. The obvious choice for this would have
 111  * been the vdwtype/coultype values in the forcerecord, but unfortunately
 112  * those types are braindead - for instance both Buckingham and normal
 113  * Lennard-Jones use the same value (evdwCUT), and a separate boolean variable
 114  * to determine which interaction is used. There is further no special value
 115  * for 'no interaction'. For backward compatibility with old TPR files we won't
 116  * change this in the 3.x series, so when calling this routine you should use:
 117  *
 118  * icoul=0 no coulomb interaction
 119  * icoul=1 cutoff standard coulomb
 120  * icoul=2 reaction-field coulomb
 121  * icoul=3 tabulated coulomb
 122  *
 123  * ivdw=0 no vdw interaction
 124  * ivdw=1 standard L-J interaction
 125  * ivdw=2 Buckingham
 126  * ivdw=3 tabulated vdw.
 127  *
 128  * Kind of ugly, but it works.
 129  */
 130 static void init_nblist(t_nblist *nl_sr,t_nblist *nl_lr,
 131                         int maxsr,int maxlr,
 132                         int ivdw, int icoul,
 133                         bool bfree, int enlist)
 134 {
 135     t_nblist *nl;
 136     int      homenr;
 137     int      i,nn;
 138
 139     int inloop[20] =
 140     {
 141         eNR_NBKERNEL_NONE,
 142         eNR_NBKERNEL010,
 143         eNR_NBKERNEL020,
 144         eNR_NBKERNEL030,
 145         eNR_NBKERNEL100,
 146         eNR_NBKERNEL110,
 147         eNR_NBKERNEL120,
 148         eNR_NBKERNEL130,
 149         eNR_NBKERNEL200,
 150         eNR_NBKERNEL210,
 151         eNR_NBKERNEL220,
 152         eNR_NBKERNEL230,
 153         eNR_NBKERNEL300,
 154         eNR_NBKERNEL310,
 155         eNR_NBKERNEL320,
 156         eNR_NBKERNEL330,
 157         eNR_NBKERNEL400,
 158         eNR_NBKERNEL410,
 159         eNR_NBKERNEL_NONE,
 160         eNR_NBKERNEL430
 161     };
 162
 163     for(i=0; (i<2); i++)
 164     {
 165         nl     = (i == 0) ? nl_sr : nl_lr;
 166         homenr = (i == 0) ? maxsr : maxlr;
 167
 168         if (nl == NULL)
 169         {
 170             continue;
 171         }
 172
 173         /* Set coul/vdw in neighborlist, and for the normal loops we determine
 174          * an index of which one to call.
 175          */
 176         nl->ivdw  = ivdw;
 177         nl->icoul = icoul;
 178         nl->free_energy = bfree;
 179
 180         if (bfree)
 181         {
 182             nl->enlist  = enlistATOM_ATOM;
 183             nl->il_code = eNR_NBKERNEL_FREE_ENERGY;
 184         }
 185         else
 186         {
 187             nl->enlist = enlist;
 188
 189             nn = inloop[4*icoul + ivdw];
 190
 191             /* solvent loops follow directly after the corresponding
 192             * ordinary loops, in the order:
 193             *
 194             * SPC, SPC-SPC, TIP4p, TIP4p-TIP4p
 195             *
 196             */
 197             switch (enlist) {
 198             case enlistATOM_ATOM:
 199             case enlistCG_CG:
 200                 break;
 201             case enlistSPC_ATOM:     nn += 1; break;
 202             case enlistSPC_SPC:      nn += 2; break;
 203             case enlistTIP4P_ATOM:   nn += 3; break;
 204             case enlistTIP4P_TIP4P:  nn += 4; break;
 205             }
 206
 207             nl->il_code = nn;
 208         }
 209
 210         if (debug)
 211             fprintf(debug,"Initiating neighbourlist type %d for %s interactions,\nwith %d SR, %d LR atoms.\n",
 212                     nl->il_code,ENLISTTYPE(enlist),maxsr,maxlr);
 213
 214         /* maxnri is influenced by the number of shifts (maximum is 8)
 215          * and the number of energy groups.
 216          * If it is not enough, nl memory will be reallocated during the run.
 217          * 4 seems to be a reasonable factor, which only causes reallocation
 218          * during runs with tiny and many energygroups.
 219          */
 220         nl->maxnri      = homenr*4;
 221         nl->maxnrj      = 0;
 222         nl->maxlen      = 0;
 223         nl->nri         = -1;
 224         nl->nrj         = 0;
 225         nl->iinr        = NULL;
 226         nl->gid         = NULL;
 227         nl->shift       = NULL;
 228         nl->jindex      = NULL;
 229         reallocate_nblist(nl);
 230         nl->jindex[0] = 0;
 231 #ifdef GMX_THREAD_SHM_FDECOMP
 232         nl->counter = 0;
 233         snew(nl->mtx,1);
 234         pthread_mutex_init(nl->mtx,NULL);
 235 #endif
 236     }
 237 }
 238
 239 void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr)
 240 {
 241    /* Make maxlr tunable! (does not seem to be a big difference though)
 242     * This parameter determines the number of i particles in a long range
 243     * neighbourlist. Too few means many function calls, too many means
 244     * cache trashing.
 245     */
 246    int maxsr,maxsr_wat,maxlr,maxlr_wat;
 247    int icoul,icoulf,ivdw;
 248    int solvent;
 249    int enlist_def,enlist_w,enlist_ww;
 250    int i;
 251    t_nblists *nbl;
 252
 253    /* maxsr     = homenr-fr->nWatMol*3; */
 254    maxsr     = homenr;
 255
 256    if (maxsr < 0)
 257    {
 258      gmx_fatal(FARGS,"%s, %d: Negative number of short range atoms.\n"
 259                  "Call your Gromacs dealer for assistance.",__FILE__,__LINE__);
 260    }
 261    /* This is just for initial allocation, so we do not reallocate
 262     * all the nlist arrays many times in a row.
 263     * The numbers seem very accurate, but they are uncritical.
 264     */
 265    maxsr_wat = min(fr->nWatMol,(homenr+2)/3);
 266    if (fr->bTwinRange)
 267    {
 268        maxlr     = 50;
 269        maxlr_wat = min(maxsr_wat,maxlr);
 270    }
 271    else
 272    {
 273      maxlr = maxlr_wat = 0;
 274    }
 275
 276    /* Determine the values for icoul/ivdw. */
 277    /* Start with GB */
 278    if(fr->bGB)
 279    {
 280        icoul=4;
 281    }
 282    else if (fr->bcoultab)
 283    {
 284        icoul = 3;
 285    }
 286    else if (EEL_RF(fr->eeltype))
 287    {
 288        icoul = 2;
 289    }
 290    else
 291    {
 292        icoul = 1;
 293    }
 294
 295    if (fr->bvdwtab)
 296    {
 297        ivdw = 3;
 298    }
 299    else if (fr->bBHAM)
 300    {
 301        ivdw = 2;
 302    }
 303    else
 304    {
 305        ivdw = 1;
 306    }
 307
 308    fr->ns.bCGlist = (getenv("GMX_NBLISTCG") != 0);
 309    if (!fr->ns.bCGlist)
 310    {
 311        enlist_def = enlistATOM_ATOM;
 312    }
 313    else
 314    {
 315        enlist_def = enlistCG_CG;
 316        if (log != NULL)
 317        {
 318            fprintf(log,"\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
 319        }
 320        if (!fr->bExcl_IntraCGAll_InterCGNone)
 321        {
 322            gmx_fatal(FARGS,"The charge-group - charge-group force loops only support systems with all intra-cg interactions excluded and no inter-cg exclusions, this is not the case for this system.");
 323        }
 324    }
 325
 326    if (fr->solvent_opt == esolTIP4P) {
 327        enlist_w  = enlistTIP4P_ATOM;
 328        enlist_ww = enlistTIP4P_TIP4P;
 329    } else {
 330        enlist_w  = enlistSPC_ATOM;
 331        enlist_ww = enlistSPC_SPC;
 332    }
 333
 334    for(i=0; i<fr->nnblists; i++)
 335    {
 336        nbl = &(fr->nblists[i]);
 337        init_nblist(&nbl->nlist_sr[eNL_VDWQQ],&nbl->nlist_lr[eNL_VDWQQ],
 338                    maxsr,maxlr,ivdw,icoul,FALSE,enlist_def);
 339        init_nblist(&nbl->nlist_sr[eNL_VDW],&nbl->nlist_lr[eNL_VDW],
 340                    maxsr,maxlr,ivdw,0,FALSE,enlist_def);
 341        init_nblist(&nbl->nlist_sr[eNL_QQ],&nbl->nlist_lr[eNL_QQ],
 342                    maxsr,maxlr,0,icoul,FALSE,enlist_def);
 343        init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATER],&nbl->nlist_lr[eNL_VDWQQ_WATER],
 344                    maxsr_wat,maxlr_wat,ivdw,icoul, FALSE,enlist_w);
 345        init_nblist(&nbl->nlist_sr[eNL_QQ_WATER],&nbl->nlist_lr[eNL_QQ_WATER],
 346                    maxsr_wat,maxlr_wat,0,icoul, FALSE,enlist_w);
 347        init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATERWATER],&nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
 348                    maxsr_wat,maxlr_wat,ivdw,icoul, FALSE,enlist_ww);
 349        init_nblist(&nbl->nlist_sr[eNL_QQ_WATERWATER],&nbl->nlist_lr[eNL_QQ_WATERWATER],
 350                    maxsr_wat,maxlr_wat,0,icoul, FALSE,enlist_ww);
 351
 352        if (fr->efep != efepNO)
 353        {
 354            if (fr->bEwald)
 355            {
 356                icoulf = 5;
 357            }
 358            else
 359            {
 360                icoulf = icoul;
 361            }
 362
 363            init_nblist(&nbl->nlist_sr[eNL_VDWQQ_FREE],&nbl->nlist_lr[eNL_VDWQQ_FREE],
 364                        maxsr,maxlr,ivdw,icoulf,TRUE,enlistATOM_ATOM);
 365            init_nblist(&nbl->nlist_sr[eNL_VDW_FREE],&nbl->nlist_lr[eNL_VDW_FREE],
 366                        maxsr,maxlr,ivdw,0,TRUE,enlistATOM_ATOM);
 367            init_nblist(&nbl->nlist_sr[eNL_QQ_FREE],&nbl->nlist_lr[eNL_QQ_FREE],
 368                        maxsr,maxlr,0,icoulf,TRUE,enlistATOM_ATOM);
 369        }
 370    }
 371    /* QMMM MM list */
 372    if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
 373    {
 374        init_nblist(&fr->QMMMlist,NULL,
 375                    maxsr,maxlr,0,icoul,FALSE,enlistATOM_ATOM);
 376    }
 377
 378    fr->ns.nblist_initialized=TRUE;
 379 }
 380
 381 static void reset_nblist(t_nblist *nl)
 382 {
 383      nl->nri       = -1;
 384      nl->nrj       = 0;
 385      nl->maxlen    = 0;
 386      if (nl->jindex)
 387      {
 388          nl->jindex[0] = 0;
 389      }
 390 }
 391
 392 static void reset_neighbor_list(t_forcerec *fr,bool bLR,int nls,int eNL)
 393 {
 394     int n,i;
 395
 396     if (bLR)
 397     {
 398         reset_nblist(&(fr->nblists[nls].nlist_lr[eNL]));
 399     }
 400     else
 401     {
 402         for(n=0; n<fr->nnblists; n++)
 403         {
 404             for(i=0; i<eNL_NR; i++)
 405             {
 406                 reset_nblist(&(fr->nblists[n].nlist_sr[i]));
 407             }
 408         }
 409         if (fr->bQMMM)
 410         {
 411             /* only reset the short-range nblist */
 412             reset_nblist(&(fr->QMMMlist));
 413         }
 414     }
 415 }
 416
 417
 418
 419
 420 static inline void new_i_nblist(t_nblist *nlist,
 421                                 bool bLR,atom_id i_atom,int shift,int gid)
 422 {
 423     int    i,k,nri,nshift;
 424
 425     nri = nlist->nri;
 426
 427     /* Check whether we have to increase the i counter */
 428     if ((nri == -1) ||
 429         (nlist->iinr[nri]  != i_atom) ||
 430         (nlist->shift[nri] != shift) ||
 431         (nlist->gid[nri]   != gid))
 432     {
 433         /* This is something else. Now see if any entries have
 434          * been added in the list of the previous atom.
 435          */
 436         if ((nri == -1) ||
 437             ((nlist->jindex[nri+1] > nlist->jindex[nri]) &&
 438              (nlist->gid[nri] != -1)))
 439         {
 440             /* If so increase the counter */
 441             nlist->nri++;
 442             nri++;
 443             if (nlist->nri >= nlist->maxnri)
 444             {
 445                 nlist->maxnri += over_alloc_large(nlist->nri);
 446                 reallocate_nblist(nlist);
 447             }
 448         }
 449         /* Set the number of neighbours and the atom number */
 450         nlist->jindex[nri+1] = nlist->jindex[nri];
 451         nlist->iinr[nri]     = i_atom;
 452         nlist->gid[nri]      = gid;
 453         nlist->shift[nri]    = shift;
 454     }
 455 }
 456
 457 static inline void close_i_nblist(t_nblist *nlist)
 458 {
 459     int nri = nlist->nri;
 460     int len;
 461
 462     if (nri >= 0)
 463     {
 464         nlist->jindex[nri+1] = nlist->nrj;
 465
 466         len=nlist->nrj -  nlist->jindex[nri];
 467
 468         /* nlist length for water i molecules is treated statically
 469          * in the innerloops
 470          */
 471         if (len > nlist->maxlen)
 472         {
 473             nlist->maxlen = len;
 474         }
 475     }
 476 }
 477
 478 static inline void close_nblist(t_nblist *nlist)
 479 {
 480     /* Only close this nblist when it has been initialized */
 481     if (nlist->jindex)
 482     {
 483         nlist->nri++;
 484     }
 485 }
 486
 487 static inline void close_neighbor_list(t_forcerec *fr,bool bLR,int nls,int eNL,
 488                                        bool bMakeQMMMnblist)
 489 {
 490     int n,i;
 491
 492     if (bMakeQMMMnblist) {
 493         if (!bLR)
 494         {
 495             close_nblist(&(fr->QMMMlist));
 496         }
 497     }
 498     else
 499     {
 500         if (bLR)
 501         {
 502             close_nblist(&(fr->nblists[nls].nlist_lr[eNL]));
 503         }
 504         else
 505         {
 506             for(n=0; n<fr->nnblists; n++)
 507             {
 508                 for(i=0; (i<eNL_NR); i++)
 509                 {
 510                     close_nblist(&(fr->nblists[n].nlist_sr[i]));
 511                 }
 512             }
 513         }
 514     }
 515 }
 516
 517 static inline void add_j_to_nblist(t_nblist *nlist,atom_id j_atom,bool bLR)
 518 {
 519     int nrj=nlist->nrj;
 520
 521     if (nlist->nrj >= nlist->maxnrj)
 522     {
 523         nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
 524         if (gmx_debug_at)
 525             fprintf(debug,"Increasing %s nblist %s j size to %d\n",
 526                     bLR ? "LR" : "SR",nrnb_str(nlist->il_code),nlist->maxnrj);
 527
 528         srenew(nlist->jjnr,nlist->maxnrj);
 529     }
 530
 531     nlist->jjnr[nrj] = j_atom;
 532     nlist->nrj ++;
 533 }
 534
 535 static inline void add_j_to_nblist_cg(t_nblist *nlist,
 536                                       atom_id j_start,int j_end,
 537                                       t_excl *bexcl,bool bLR)
 538 {
 539     int nrj=nlist->nrj;
 540     int j;
 541
 542     if (nlist->nrj >= nlist->maxnrj)
 543     {
 544         nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
 545         if (gmx_debug_at)
 546             fprintf(debug,"Increasing %s nblist %s j size to %d\n",
 547                     bLR ? "LR" : "SR",nrnb_str(nlist->il_code),nlist->maxnrj);
 548
 549         srenew(nlist->jjnr    ,nlist->maxnrj);
 550         srenew(nlist->jjnr_end,nlist->maxnrj);
 551         srenew(nlist->excl    ,nlist->maxnrj*MAX_CGCGSIZE);
 552     }
 553
 554     nlist->jjnr[nrj]     = j_start;
 555     nlist->jjnr_end[nrj] = j_end;
 556
 557     if (j_end - j_start > MAX_CGCGSIZE)
 558     {
 559         gmx_fatal(FARGS,"The charge-group - charge-group neighborlist do not support charge groups larger than %d, found a charge group of size %d",MAX_CGCGSIZE,j_end-j_start);
 560     }
 561
 562     /* Set the exclusions */
 563     for(j=j_start; j<j_end; j++)
 564     {
 565         nlist->excl[nrj*MAX_CGCGSIZE + j - j_start] = bexcl[j];
 566     }
 567
 568     nlist->nrj ++;
 569 }
 570
 571 typedef void
 572 put_in_list_t(bool              bHaveVdW[],
 573               int               ngid,
 574               t_mdatoms *       md,
 575               int               icg,
 576               int               jgid,
 577               int               nj,
 578               atom_id           jjcg[],
 579               atom_id           index[],
 580               t_excl            bExcl[],
 581               int               shift,
 582               t_forcerec *      fr,
 583               bool              bLR,
 584               bool              bDoVdW,
 585               bool              bDoCoul);
 586
 587 static void
 588 put_in_list_at(bool              bHaveVdW[],
 589                int               ngid,
 590                t_mdatoms *       md,
 591                int               icg,
 592                int               jgid,
 593                int               nj,
 594                atom_id           jjcg[],
 595                atom_id           index[],
 596                t_excl            bExcl[],
 597                int               shift,
 598                t_forcerec *      fr,
 599                bool              bLR,
 600                bool              bDoVdW,
 601                bool              bDoCoul)
 602 {
 603     /* The a[] index has been removed,
 604      * to put it back in i_atom should be a[i0] and jj should be a[jj].
 605      */
 606     t_nblist *   vdwc;
 607     t_nblist *   vdw;
 608     t_nblist *   coul;
 609     t_nblist *   vdwc_free  = NULL;
 610     t_nblist *   vdw_free   = NULL;
 611     t_nblist *   coul_free  = NULL;
 612     t_nblist *   vdwc_ww    = NULL;
 613     t_nblist *   coul_ww    = NULL;
 614
 615     int             i,j,jcg,igid,gid,nbl_ind,ind_ij;
 616     atom_id   jj,jj0,jj1,i_atom;
 617     int       i0,nicg,len;
 618
 619     int       *cginfo;
 620     int       *type,*typeB;
 621     real      *charge,*chargeB;
 622     real      qi,qiB,qq,rlj;
 623     bool      bFreeEnergy,bFree,bFreeJ,bNotEx,*bPert;
 624     bool      bDoVdW_i,bDoCoul_i,bDoCoul_i_sol;
 625     int       iwater,jwater;
 626     t_nblist  *nlist;
 627
 628     /* Copy some pointers */
 629     cginfo  = fr->cginfo;
 630     charge  = md->chargeA;
 631     chargeB = md->chargeB;
 632     type    = md->typeA;
 633     typeB   = md->typeB;
 634     bPert   = md->bPerturbed;
 635
 636     /* Get atom range */
 637     i0     = index[icg];
 638     nicg   = index[icg+1]-i0;
 639
 640     /* Get the i charge group info */
 641     igid   = GET_CGINFO_GID(cginfo[icg]);
 642     iwater = GET_CGINFO_SOLOPT(cginfo[icg]);
 643
 644     bFreeEnergy = FALSE;
 645     if (md->nPerturbed)
 646     {
 647         /* Check if any of the particles involved are perturbed.
 648          * If not we can do the cheaper normal put_in_list
 649          * and use more solvent optimization.
 650          */
 651         for(i=0; i<nicg; i++)
 652         {
 653             bFreeEnergy |= bPert[i0+i];
 654         }
 655         /* Loop over the j charge groups */
 656         for(j=0; (j<nj && !bFreeEnergy); j++)
 657         {
 658             jcg = jjcg[j];
 659             jj0 = index[jcg];
 660             jj1 = index[jcg+1];
 661             /* Finally loop over the atoms in the j-charge group */
 662             for(jj=jj0; jj<jj1; jj++)
 663             {
 664                 bFreeEnergy |= bPert[jj];
 665             }
 666         }
 667     }
 668
 669     /* Unpack pointers to neighbourlist structs */
 670     if (fr->nnblists == 1)
 671     {
 672         nbl_ind = 0;
 673     }
 674     else
 675     {
 676         nbl_ind = fr->gid2nblists[GID(igid,jgid,ngid)];
 677     }
 678     if (bLR)
 679     {
 680         nlist = fr->nblists[nbl_ind].nlist_lr;
 681     }
 682     else
 683     {
 684         nlist = fr->nblists[nbl_ind].nlist_sr;
 685     }
 686
 687     if (iwater != esolNO)
 688     {
 689         vdwc = &nlist[eNL_VDWQQ_WATER];
 690         vdw  = &nlist[eNL_VDW];
 691         coul = &nlist[eNL_QQ_WATER];
 692 #ifndef DISABLE_WATERWATER_NLIST
 693         vdwc_ww = &nlist[eNL_VDWQQ_WATERWATER];
 694         coul_ww = &nlist[eNL_QQ_WATERWATER];
 695 #endif
 696     }
 697     else
 698     {
 699         vdwc = &nlist[eNL_VDWQQ];
 700         vdw  = &nlist[eNL_VDW];
 701         coul = &nlist[eNL_QQ];
 702     }
 703
 704     if (!bFreeEnergy)
 705     {
 706         if (iwater != esolNO)
 707         {
 708             /* Loop over the atoms in the i charge group */
 709             i_atom  = i0;
 710             gid     = GID(igid,jgid,ngid);
 711             /* Create new i_atom for each energy group */
 712             if (bDoCoul && bDoVdW)
 713             {
 714                 new_i_nblist(vdwc,bLR,i_atom,shift,gid);
 715 #ifndef DISABLE_WATERWATER_NLIST
 716                 new_i_nblist(vdwc_ww,bLR,i_atom,shift,gid);
 717 #endif
 718             }
 719             if (bDoVdW)
 720             {
 721                 new_i_nblist(vdw,bLR,i_atom,shift,gid);
 722             }
 723             if (bDoCoul)
 724             {
 725                 new_i_nblist(coul,bLR,i_atom,shift,gid);
 726 #ifndef DISABLE_WATERWATER_NLIST
 727                 new_i_nblist(coul_ww,bLR,i_atom,shift,gid);
 728 #endif
 729             }
 730           /* Loop over the j charge groups */
 731             for(j=0; (j<nj); j++)
 732             {
 733                 jcg=jjcg[j];
 734
 735                 if (jcg == icg)
 736                 {
 737                     continue;
 738                 }
 739
 740                 jj0 = index[jcg];
 741                 jwater = GET_CGINFO_SOLOPT(cginfo[jcg]);
 742
 743                 if (iwater == esolSPC && jwater == esolSPC)
 744                 {
 745                     /* Interaction between two SPC molecules */
 746                     if (!bDoCoul)
 747                     {
 748                         /* VdW only - only first atoms in each water interact */
 749                         add_j_to_nblist(vdw,jj0,bLR);
 750                     }
 751                     else
 752                     {
 753 #ifdef DISABLE_WATERWATER_NLIST
 754                         /* Add entries for the three atoms - only do VdW if we need to */
 755                         if (!bDoVdW)
 756                         {
 757                             add_j_to_nblist(coul,jj0,bLR);
 758                         }
 759                         else
 760                         {
 761                             add_j_to_nblist(vdwc,jj0,bLR);
 762                         }
 763                         add_j_to_nblist(coul,jj0+1,bLR);
 764                         add_j_to_nblist(coul,jj0+2,bLR);
 765 #else
 766                         /* One entry for the entire water-water interaction */
 767                         if (!bDoVdW)
 768                         {
 769                             add_j_to_nblist(coul_ww,jj0,bLR);
 770                         }
 771                         else
 772                         {
 773                             add_j_to_nblist(vdwc_ww,jj0,bLR);
 774                         }
 775 #endif
 776                     }
 777                 }
 778                 else if (iwater == esolTIP4P && jwater == esolTIP4P)
 779                 {
 780                     /* Interaction between two TIP4p molecules */
 781                     if (!bDoCoul)
 782                     {
 783                         /* VdW only - only first atoms in each water interact */
 784                         add_j_to_nblist(vdw,jj0,bLR);
 785                     }
 786                     else
 787                     {
 788 #ifdef DISABLE_WATERWATER_NLIST
 789                         /* Add entries for the four atoms - only do VdW if we need to */
 790                         if (bDoVdW)
 791                         {
 792                             add_j_to_nblist(vdw,jj0,bLR);
 793                         }
 794                         add_j_to_nblist(coul,jj0+1,bLR);
 795                         add_j_to_nblist(coul,jj0+2,bLR);
 796                         add_j_to_nblist(coul,jj0+3,bLR);
 797 #else
 798                         /* One entry for the entire water-water interaction */
 799                         if (!bDoVdW)
 800                         {
 801                             add_j_to_nblist(coul_ww,jj0,bLR);
 802                         }
 803                         else
 804                         {
 805                             add_j_to_nblist(vdwc_ww,jj0,bLR);
 806                         }
 807 #endif
 808                     }
 809                 }
 810                 else
 811                 {
 812                     /* j charge group is not water, but i is.
 813                      * Add entries to the water-other_atom lists; the geometry of the water
 814                      * molecule doesn't matter - that is taken care of in the nonbonded kernel,
 815                      * so we don't care if it is SPC or TIP4P...
 816                      */
 817
 818                     jj1 = index[jcg+1];
 819
 820                     if (!bDoVdW)
 821                     {
 822                         for(jj=jj0; (jj<jj1); jj++)
 823                         {
 824                             if (charge[jj] != 0)
 825                             {
 826                                 add_j_to_nblist(coul,jj,bLR);
 827                             }
 828                         }
 829                     }
 830                     else if (!bDoCoul)
 831                     {
 832                         for(jj=jj0; (jj<jj1); jj++)
 833                         {
 834                             if (bHaveVdW[type[jj]])
 835                             {
 836                                 add_j_to_nblist(vdw,jj,bLR);
 837                             }
 838                         }
 839                     }
 840                     else
 841                     {
 842                         /* _charge_ _groups_ interact with both coulomb and LJ */
 843                         /* Check which atoms we should add to the lists!       */
 844                         for(jj=jj0; (jj<jj1); jj++)
 845                         {
 846                             if (bHaveVdW[type[jj]])
 847                             {
 848                                 if (charge[jj] != 0)
 849                                 {
 850                                     add_j_to_nblist(vdwc,jj,bLR);
 851                                 }
 852                                 else
 853                                 {
 854                                     add_j_to_nblist(vdw,jj,bLR);
 855                                 }
 856                             }
 857                             else if (charge[jj] != 0)
 858                             {
 859                                 add_j_to_nblist(coul,jj,bLR);
 860                             }
 861                         }
 862                     }
 863                 }
 864             }
 865             close_i_nblist(vdw);
 866             close_i_nblist(coul);
 867             close_i_nblist(vdwc);
 868 #ifndef DISABLE_WATERWATER_NLIST
 869             close_i_nblist(coul_ww);
 870             close_i_nblist(vdwc_ww);
 871 #endif
 872         }
 873         else
 874         {
 875             /* no solvent as i charge group */
 876             /* Loop over the atoms in the i charge group */
 877             for(i=0; i<nicg; i++)
 878             {
 879                 i_atom  = i0+i;
 880                 gid     = GID(igid,jgid,ngid);
 881                 qi      = charge[i_atom];
 882
 883                 /* Create new i_atom for each energy group */
 884                 if (bDoVdW && bDoCoul)
 885                 {
 886                     new_i_nblist(vdwc,bLR,i_atom,shift,gid);
 887                 }
 888                 if (bDoVdW)
 889                 {
 890                     new_i_nblist(vdw,bLR,i_atom,shift,gid);
 891                 }
 892                 if (bDoCoul)
 893                 {
 894                     new_i_nblist(coul,bLR,i_atom,shift,gid);
 895                 }
 896                 bDoVdW_i  = (bDoVdW  && bHaveVdW[type[i_atom]]);
 897                 bDoCoul_i = (bDoCoul && qi!=0);
 898
 899                 if (bDoVdW_i || bDoCoul_i)
 900                 {
 901                     /* Loop over the j charge groups */
 902                     for(j=0; (j<nj); j++)
 903                     {
 904                         jcg=jjcg[j];
 905
 906                         /* Check for large charge groups */
 907                         if (jcg == icg)
 908                         {
 909                             jj0 = i0 + i + 1;
 910                         }
 911                         else
 912                         {
 913                             jj0 = index[jcg];
 914                         }
 915
 916                         jj1=index[jcg+1];
 917                         /* Finally loop over the atoms in the j-charge group */
 918                         for(jj=jj0; jj<jj1; jj++)
 919                         {
 920                             bNotEx = NOTEXCL(bExcl,i,jj);
 921
 922                             if (bNotEx)
 923                             {
 924                                 if (!bDoVdW_i)
 925                                 {
 926                                     if (charge[jj] != 0)
 927                                     {
 928                                         add_j_to_nblist(coul,jj,bLR);
 929                                     }
 930                                 }
 931                                 else if (!bDoCoul_i)
 932                                 {
 933                                     if (bHaveVdW[type[jj]])
 934                                     {
 935                                         add_j_to_nblist(vdw,jj,bLR);
 936                                     }
 937                                 }
 938                                 else
 939                                 {
 940                                     if (bHaveVdW[type[jj]])
 941                                     {
 942                                         if (charge[jj] != 0)
 943                                         {
 944                                             add_j_to_nblist(vdwc,jj,bLR);
 945                                         }
 946                                         else
 947                                         {
 948                                             add_j_to_nblist(vdw,jj,bLR);
 949                                         }
 950                                     }
 951                                     else if (charge[jj] != 0)
 952                                     {
 953                                         add_j_to_nblist(coul,jj,bLR);
 954                                     }
 955                                 }
 956                             }
 957                         }
 958                     }
 959                 }
 960                 close_i_nblist(vdw);
 961                 close_i_nblist(coul);
 962                 close_i_nblist(vdwc);
 963             }
 964         }
 965     }
 966     else
 967     {
 968         /* we are doing free energy */
 969         vdwc_free = &nlist[eNL_VDWQQ_FREE];
 970         vdw_free  = &nlist[eNL_VDW_FREE];
 971         coul_free = &nlist[eNL_QQ_FREE];
 972         /* Loop over the atoms in the i charge group */
 973         for(i=0; i<nicg; i++)
 974         {
 975             i_atom  = i0+i;
 976             gid     = GID(igid,jgid,ngid);
 977             qi      = charge[i_atom];
 978             qiB     = chargeB[i_atom];
 979
 980             /* Create new i_atom for each energy group */
 981             if (bDoVdW && bDoCoul)
 982                 new_i_nblist(vdwc,bLR,i_atom,shift,gid);
 983             if (bDoVdW)
 984                 new_i_nblist(vdw,bLR,i_atom,shift,gid);
 985             if (bDoCoul)
 986                 new_i_nblist(coul,bLR,i_atom,shift,gid);
 987
 988             new_i_nblist(vdw_free,bLR,i_atom,shift,gid);
 989             new_i_nblist(coul_free,bLR,i_atom,shift,gid);
 990             new_i_nblist(vdwc_free,bLR,i_atom,shift,gid);
 991
 992             bDoVdW_i  = (bDoVdW  &&
 993                          (bHaveVdW[type[i_atom]] || bHaveVdW[typeB[i_atom]]));
 994             bDoCoul_i = (bDoCoul && (qi!=0 || qiB!=0));
 995             /* For TIP4P the first atom does not have a charge,
 996              * but the last three do. So we should still put an atom
 997              * without LJ but with charge in the water-atom neighborlist
 998              * for a TIP4p i charge group.
 999              * For SPC type water the first atom has LJ and charge,
1000              * so there is no such problem.
1001              */
1002             if (iwater == esolNO)
1003             {
1004                 bDoCoul_i_sol = bDoCoul_i;
1005             }
1006             else
1007             {
1008                 bDoCoul_i_sol = bDoCoul;
1009             }
1010
1011             if (bDoVdW_i || bDoCoul_i_sol)
1012             {
1013                 /* Loop over the j charge groups */
1014                 for(j=0; (j<nj); j++)
1015                 {
1016                     jcg=jjcg[j];
1017
1018                     /* Check for large charge groups */
1019                     if (jcg == icg)
1020                     {
1021                         jj0 = i0 + i + 1;
1022                     }
1023                     else
1024                     {
1025                         jj0 = index[jcg];
1026                     }
1027
1028                     jj1=index[jcg+1];
1029                     /* Finally loop over the atoms in the j-charge group */
1030                     bFree = bPert[i_atom];
1031                     for(jj=jj0; (jj<jj1); jj++)
1032                     {
1033                         bFreeJ = bFree || bPert[jj];
1034                         /* Complicated if, because the water H's should also
1035                          * see perturbed j-particles
1036                          */
1037                         if (iwater==esolNO || i==0 || bFreeJ)
1038                         {
1039                             bNotEx = NOTEXCL(bExcl,i,jj);
1040
1041                             if (bNotEx)
1042                             {
1043                                 if (bFreeJ)
1044                                 {
1045                                     if (!bDoVdW_i)
1046                                     {
1047                                         if (charge[jj]!=0 || chargeB[jj]!=0)
1048                                         {
1049                                             add_j_to_nblist(coul_free,jj,bLR);
1050                                         }
1051                                     }
1052                                     else if (!bDoCoul_i)
1053                                     {
1054                                         if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
1055                                         {
1056                                             add_j_to_nblist(vdw_free,jj,bLR);
1057                                         }
1058                                     }
1059                                     else
1060                                     {
1061                                         if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
1062                                         {
1063                                             if (charge[jj]!=0 || chargeB[jj]!=0)
1064                                             {
1065                                                 add_j_to_nblist(vdwc_free,jj,bLR);
1066                                             }
1067                                             else
1068                                             {
1069                                                 add_j_to_nblist(vdw_free,jj,bLR);
1070                                             }
1071                                         }
1072                                         else if (charge[jj]!=0 || chargeB[jj]!=0)
1073                                             add_j_to_nblist(coul_free,jj,bLR);
1074                                     }
1075                                 }
1076                                 else if (!bDoVdW_i)
1077                                 {
1078                                     /* This is done whether or not bWater is set */
1079                                     if (charge[jj] != 0)
1080                                     {
1081                                         add_j_to_nblist(coul,jj,bLR);
1082                                     }
1083                                 }
1084                                 else if (!bDoCoul_i_sol)
1085                                 {
1086                                     if (bHaveVdW[type[jj]])
1087                                     {
1088                                         add_j_to_nblist(vdw,jj,bLR);
1089                                     }
1090                                 }
1091                                 else
1092                                 {
1093                                     if (bHaveVdW[type[jj]])
1094                                     {
1095                                         if (charge[jj] != 0)
1096                                         {
1097                                             add_j_to_nblist(vdwc,jj,bLR);
1098                                         }
1099                                         else
1100                                         {
1101                                             add_j_to_nblist(vdw,jj,bLR);
1102                                         }
1103                                     }
1104                                     else if (charge[jj] != 0)
1105                                     {
1106                                         add_j_to_nblist(coul,jj,bLR);
1107                                     }
1108                                 }
1109                             }
1110                         }
1111                     }
1112                 }
1113             }
1114             close_i_nblist(vdw);
1115             close_i_nblist(coul);
1116             close_i_nblist(vdwc);
1117             close_i_nblist(vdw_free);
1118             close_i_nblist(coul_free);
1119             close_i_nblist(vdwc_free);
1120         }
1121     }
1122 }
1123
1124 static void
1125 put_in_list_qmmm(bool              bHaveVdW[],
1126                  int               ngid,
1127                  t_mdatoms *       md,
1128                  int               icg,
1129                  int               jgid,
1130                  int               nj,
1131                  atom_id           jjcg[],
1132                  atom_id           index[],
1133                  t_excl            bExcl[],
1134                  int               shift,
1135                  t_forcerec *      fr,
1136                  bool              bLR,
1137                  bool              bDoVdW,
1138                  bool              bDoCoul)
1139 {
1140     t_nblist *   coul;
1141     int           i,j,jcg,igid,gid;
1142     atom_id   jj,jj0,jj1,i_atom;
1143     int       i0,nicg;
1144     bool      bNotEx;
1145
1146     /* Get atom range */
1147     i0     = index[icg];
1148     nicg   = index[icg+1]-i0;
1149
1150     /* Get the i charge group info */
1151     igid   = GET_CGINFO_GID(fr->cginfo[icg]);
1152
1153     coul = &fr->QMMMlist;
1154
1155     /* Loop over atoms in the ith charge group */
1156     for (i=0;i<nicg;i++)
1157     {
1158         i_atom = i0+i;
1159         gid    = GID(igid,jgid,ngid);
1160         /* Create new i_atom for each energy group */
1161         new_i_nblist(coul,bLR,i_atom,shift,gid);
1162
1163         /* Loop over the j charge groups */
1164         for (j=0;j<nj;j++)
1165         {
1166             jcg=jjcg[j];
1167
1168             /* Charge groups cannot have QM and MM atoms simultaneously */
1169             if (jcg!=icg)
1170             {
1171                 jj0 = index[jcg];
1172                 jj1 = index[jcg+1];
1173                 /* Finally loop over the atoms in the j-charge group */
1174                 for(jj=jj0; jj<jj1; jj++)
1175                 {
1176                     bNotEx = NOTEXCL(bExcl,i,jj);
1177                     if(bNotEx)
1178                         add_j_to_nblist(coul,jj,bLR);
1179                 }
1180             }
1181         }
1182         close_i_nblist(coul);
1183     }
1184 }
1185
1186 static void
1187 put_in_list_cg(bool              bHaveVdW[],
1188                int               ngid,
1189                t_mdatoms *       md,
1190                int               icg,
1191                int               jgid,
1192                int               nj,
1193                atom_id           jjcg[],
1194                atom_id           index[],
1195                t_excl            bExcl[],
1196                int               shift,
1197                t_forcerec *      fr,
1198                bool              bLR,
1199                bool              bDoVdW,
1200                bool              bDoCoul)
1201 {
1202     int          cginfo;
1203     int          igid,gid,nbl_ind;
1204     t_nblist *   vdwc;
1205     int          j,jcg;
1206
1207     cginfo = fr->cginfo[icg];
1208
1209     igid = GET_CGINFO_GID(cginfo);
1210     gid  = GID(igid,jgid,ngid);
1211
1212     /* Unpack pointers to neighbourlist structs */
1213     if (fr->nnblists == 1)
1214     {
1215         nbl_ind = 0;
1216     }
1217     else
1218     {
1219         nbl_ind = fr->gid2nblists[gid];
1220     }
1221     if (bLR)
1222     {
1223         vdwc = &fr->nblists[nbl_ind].nlist_lr[eNL_VDWQQ];
1224     }
1225     else
1226     {
1227         vdwc = &fr->nblists[nbl_ind].nlist_sr[eNL_VDWQQ];
1228     }
1229
1230     /* Make a new neighbor list for charge group icg.
1231      * Currently simply one neighbor list is made with LJ and Coulomb.
1232      * If required, zero interactions could be removed here
1233      * or in the force loop.
1234      */
1235     new_i_nblist(vdwc,bLR,index[icg],shift,gid);
1236     vdwc->iinr_end[vdwc->nri] = index[icg+1];
1237
1238     for(j=0; (j<nj); j++)
1239     {
1240         jcg = jjcg[j];
1241         /* Skip the icg-icg pairs if all self interactions are excluded */
1242         if (!(jcg == icg && GET_CGINFO_EXCL_INTRA(cginfo)))
1243         {
1244             /* Here we add the j charge group jcg to the list,
1245              * exclusions are also added to the list.
1246              */
1247             add_j_to_nblist_cg(vdwc,index[jcg],index[jcg+1],bExcl,bLR);
1248         }
1249     }
1250
1251     close_i_nblist(vdwc);
1252 }
1253
1254 static void setexcl(atom_id start,atom_id end,t_blocka *excl,bool b,
1255                     t_excl bexcl[])
1256 {
1257     atom_id i,k;
1258
1259     if (b)
1260     {
1261         for(i=start; i<end; i++)
1262         {
1263             for(k=excl->index[i]; k<excl->index[i+1]; k++)
1264             {
1265                 SETEXCL(bexcl,i-start,excl->a[k]);
1266             }
1267         }
1268     }
1269     else
1270     {
1271         for(i=start; i<end; i++)
1272         {
1273             for(k=excl->index[i]; k<excl->index[i+1]; k++)
1274             {
1275                 RMEXCL(bexcl,i-start,excl->a[k]);
1276             }
1277         }
1278     }
1279 }
1280
1281 int calc_naaj(int icg,int cgtot)
1282 {
1283     int naaj;
1284
1285     if ((cgtot % 2) == 1)
1286     {
1287         /* Odd number of charge groups, easy */
1288         naaj = 1 + (cgtot/2);
1289     }
1290     else if ((cgtot % 4) == 0)
1291     {
1292     /* Multiple of four is hard */
1293         if (icg < cgtot/2)
1294         {
1295             if ((icg % 2) == 0)
1296             {
1297                 naaj=1+(cgtot/2);
1298             }
1299             else
1300             {
1301                 naaj=cgtot/2;
1302             }
1303         }
1304         else
1305         {
1306             if ((icg % 2) == 1)
1307             {
1308                 naaj=1+(cgtot/2);
1309             }
1310             else
1311             {
1312                 naaj=cgtot/2;
1313             }
1314         }
1315     }
1316     else
1317     {
1318         /* cgtot/2 = odd */
1319         if ((icg % 2) == 0)
1320         {
1321             naaj=1+(cgtot/2);
1322         }
1323         else
1324         {
1325             naaj=cgtot/2;
1326         }
1327     }
1328 #ifdef DEBUG
1329     fprintf(log,"naaj=%d\n",naaj);
1330 #endif
1331
1332     return naaj;
1333 }
1334
1335 /************************************************
1336  *
1337  *  S I M P L E      C O R E     S T U F F
1338  *
1339  ************************************************/
1340
1341 static real calc_image_tric(rvec xi,rvec xj,matrix box,
1342                             rvec b_inv,int *shift)
1343 {
1344     /* This code assumes that the cut-off is smaller than
1345      * a half times the smallest diagonal element of the box.
1346      */
1347     const real h25=2.5;
1348     real dx,dy,dz;
1349     real r2;
1350     int  tx,ty,tz;
1351
1352     /* Compute diff vector */
1353     dz = xj[ZZ] - xi[ZZ];
1354     dy = xj[YY] - xi[YY];
1355     dx = xj[XX] - xi[XX];
1356
1357   /* Perform NINT operation, using trunc operation, therefore
1358    * we first add 2.5 then subtract 2 again
1359    */
1360     tz = dz*b_inv[ZZ] + h25;
1361     tz -= 2;
1362     dz -= tz*box[ZZ][ZZ];
1363     dy -= tz*box[ZZ][YY];
1364     dx -= tz*box[ZZ][XX];
1365
1366     ty = dy*b_inv[YY] + h25;
1367     ty -= 2;
1368     dy -= ty*box[YY][YY];
1369     dx -= ty*box[YY][XX];
1370
1371     tx = dx*b_inv[XX]+h25;
1372     tx -= 2;
1373     dx -= tx*box[XX][XX];
1374
1375     /* Distance squared */
1376     r2 = (dx*dx) + (dy*dy) + (dz*dz);
1377
1378     *shift = XYZ2IS(tx,ty,tz);
1379
1380     return r2;
1381 }
1382
1383 static real calc_image_rect(rvec xi,rvec xj,rvec box_size,
1384                             rvec b_inv,int *shift)
1385 {
1386     const real h15=1.5;
1387     real ddx,ddy,ddz;
1388     real dx,dy,dz;
1389     real r2;
1390     int  tx,ty,tz;
1391
1392     /* Compute diff vector */
1393     dx = xj[XX] - xi[XX];
1394     dy = xj[YY] - xi[YY];
1395     dz = xj[ZZ] - xi[ZZ];
1396
1397     /* Perform NINT operation, using trunc operation, therefore
1398      * we first add 1.5 then subtract 1 again
1399      */
1400     tx = dx*b_inv[XX] + h15;
1401     ty = dy*b_inv[YY] + h15;
1402     tz = dz*b_inv[ZZ] + h15;
1403     tx--;
1404     ty--;
1405     tz--;
1406
1407     /* Correct diff vector for translation */
1408     ddx = tx*box_size[XX] - dx;
1409     ddy = ty*box_size[YY] - dy;
1410     ddz = tz*box_size[ZZ] - dz;
1411
1412     /* Distance squared */
1413     r2 = (ddx*ddx) + (ddy*ddy) + (ddz*ddz);
1414
1415     *shift = XYZ2IS(tx,ty,tz);
1416
1417     return r2;
1418 }
1419
1420 static void add_simple(t_ns_buf *nsbuf,int nrj,atom_id cg_j,
1421                        bool bHaveVdW[],int ngid,t_mdatoms *md,
1422                        int icg,int jgid,t_block *cgs,t_excl bexcl[],
1423                        int shift,t_forcerec *fr,put_in_list_t *put_in_list)
1424 {
1425     if (nsbuf->nj + nrj > MAX_CG)
1426     {
1427         put_in_list(bHaveVdW,ngid,md,icg,jgid,nsbuf->ncg,nsbuf->jcg,
1428                     cgs->index,bexcl,shift,fr,FALSE,TRUE,TRUE);
1429         /* Reset buffer contents */
1430         nsbuf->ncg = nsbuf->nj = 0;
1431     }
1432     nsbuf->jcg[nsbuf->ncg++] = cg_j;
1433     nsbuf->nj += nrj;
1434 }
1435
1436 static void ns_inner_tric(rvec x[],int icg,int *i_egp_flags,
1437                           int njcg,atom_id jcg[],
1438                           matrix box,rvec b_inv,real rcut2,
1439                           t_block *cgs,t_ns_buf **ns_buf,
1440                           bool bHaveVdW[],int ngid,t_mdatoms *md,
1441                           t_excl bexcl[],t_forcerec *fr,
1442                           put_in_list_t *put_in_list)
1443 {
1444     int      shift;
1445     int      j,nrj,jgid;
1446     int      *cginfo=fr->cginfo;
1447     atom_id  cg_j,*cgindex;
1448     t_ns_buf *nsbuf;
1449
1450     cgindex = cgs->index;
1451     shift   = CENTRAL;
1452     for(j=0; (j<njcg); j++)
1453     {
1454         cg_j   = jcg[j];
1455         nrj    = cgindex[cg_j+1]-cgindex[cg_j];
1456         if (calc_image_tric(x[icg],x[cg_j],box,b_inv,&shift) < rcut2)
1457         {
1458             jgid  = GET_CGINFO_GID(cginfo[cg_j]);
1459             if (!(i_egp_flags[jgid] & EGP_EXCL))
1460             {
1461                 add_simple(&ns_buf[jgid][shift],nrj,cg_j,
1462                            bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,shift,fr,
1463                            put_in_list);
1464             }
1465         }
1466     }
1467 }
1468
1469 static void ns_inner_rect(rvec x[],int icg,int *i_egp_flags,
1470                           int njcg,atom_id jcg[],
1471                           bool bBox,rvec box_size,rvec b_inv,real rcut2,
1472                           t_block *cgs,t_ns_buf **ns_buf,
1473                           bool bHaveVdW[],int ngid,t_mdatoms *md,
1474                           t_excl bexcl[],t_forcerec *fr,
1475                           put_in_list_t *put_in_list)
1476 {
1477     int      shift;
1478     int      j,nrj,jgid;
1479     int      *cginfo=fr->cginfo;
1480     atom_id  cg_j,*cgindex;
1481     t_ns_buf *nsbuf;
1482
1483     cgindex = cgs->index;
1484     if (bBox)
1485     {
1486         shift = CENTRAL;
1487         for(j=0; (j<njcg); j++)
1488         {
1489             cg_j   = jcg[j];
1490             nrj    = cgindex[cg_j+1]-cgindex[cg_j];
1491             if (calc_image_rect(x[icg],x[cg_j],box_size,b_inv,&shift) < rcut2)
1492             {
1493                 jgid  = GET_CGINFO_GID(cginfo[cg_j]);
1494                 if (!(i_egp_flags[jgid] & EGP_EXCL))
1495                 {
1496                     add_simple(&ns_buf[jgid][shift],nrj,cg_j,
1497                                bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,shift,fr,
1498                                put_in_list);
1499                 }
1500             }
1501         }
1502     }
1503     else
1504     {
1505         for(j=0; (j<njcg); j++)
1506         {
1507             cg_j   = jcg[j];
1508             nrj    = cgindex[cg_j+1]-cgindex[cg_j];
1509             if ((rcut2 == 0) || (distance2(x[icg],x[cg_j]) < rcut2)) {
1510                 jgid  = GET_CGINFO_GID(cginfo[cg_j]);
1511                 if (!(i_egp_flags[jgid] & EGP_EXCL))
1512                 {
1513                     add_simple(&ns_buf[jgid][CENTRAL],nrj,cg_j,
1514                                bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,CENTRAL,fr,
1515                                put_in_list);
1516                 }
1517             }
1518         }
1519     }
1520 }
1521
1522 /* ns_simple_core needs to be adapted for QMMM still 2005 */
1523
1524 static int ns_simple_core(t_forcerec *fr,
1525                           gmx_localtop_t *top,
1526                           t_mdatoms *md,
1527                           matrix box,rvec box_size,
1528                           t_excl bexcl[],atom_id *aaj,
1529                           int ngid,t_ns_buf **ns_buf,
1530                           put_in_list_t *put_in_list,bool bHaveVdW[])
1531 {
1532     int      naaj,k;
1533     real     rlist2;
1534     int      nsearch,icg,jcg,igid,i0,nri,nn;
1535     int      *cginfo;
1536     t_ns_buf *nsbuf;
1537     /* atom_id  *i_atoms; */
1538     t_block  *cgs=&(top->cgs);
1539     t_blocka *excl=&(top->excls);
1540     rvec     b_inv;
1541     int      m;
1542     bool     bBox,bTriclinic;
1543     int      *i_egp_flags;
1544
1545     rlist2 = sqr(fr->rlist);
1546
1547     bBox = (fr->ePBC != epbcNONE);
1548     if (bBox)
1549     {
1550         for(m=0; (m<DIM); m++)
1551         {
1552             b_inv[m] = divide(1.0,box_size[m]);
1553         }
1554         bTriclinic = TRICLINIC(box);
1555     }
1556     else
1557     {
1558         bTriclinic = FALSE;
1559     }
1560
1561     cginfo = fr->cginfo;
1562
1563     nsearch=0;
1564     for (icg=fr->cg0; (icg<fr->hcg); icg++)
1565     {
1566         /*
1567           i0        = cgs->index[icg];
1568           nri       = cgs->index[icg+1]-i0;
1569           i_atoms   = &(cgs->a[i0]);
1570           i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
1571           setexcl(nri,i_atoms,excl,TRUE,bexcl);
1572         */
1573         igid = GET_CGINFO_GID(cginfo[icg]);
1574         i_egp_flags = fr->egp_flags + ngid*igid;
1575         setexcl(cgs->index[icg],cgs->index[icg+1],excl,TRUE,bexcl);
1576
1577         naaj=calc_naaj(icg,cgs->nr);
1578         if (bTriclinic)
1579         {
1580             ns_inner_tric(fr->cg_cm,icg,i_egp_flags,naaj,&(aaj[icg]),
1581                           box,b_inv,rlist2,cgs,ns_buf,
1582                           bHaveVdW,ngid,md,bexcl,fr,put_in_list);
1583         }
1584         else
1585         {
1586             ns_inner_rect(fr->cg_cm,icg,i_egp_flags,naaj,&(aaj[icg]),
1587                           bBox,box_size,b_inv,rlist2,cgs,ns_buf,
1588                           bHaveVdW,ngid,md,bexcl,fr,put_in_list);
1589         }
1590         nsearch += naaj;
1591
1592         for(nn=0; (nn<ngid); nn++)
1593         {
1594             for(k=0; (k<SHIFTS); k++)
1595             {
1596                 nsbuf = &(ns_buf[nn][k]);
1597                 if (nsbuf->ncg > 0)
1598                 {
1599                     put_in_list(bHaveVdW,ngid,md,icg,nn,nsbuf->ncg,nsbuf->jcg,
1600                                 cgs->index,bexcl,k,fr,FALSE,TRUE,TRUE);
1601                     nsbuf->ncg=nsbuf->nj=0;
1602                 }
1603             }
1604         }
1605         /* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
1606         setexcl(cgs->index[icg],cgs->index[icg+1],excl,FALSE,bexcl);
1607     }
1608     close_neighbor_list(fr,FALSE,-1,-1,FALSE);
1609
1610     return nsearch;
1611 }
1612
1613 /************************************************
1614  *
1615  *    N S 5     G R I D     S T U F F
1616  *
1617  ************************************************/
1618
1619 static inline void get_dx(int Nx,real gridx,real rc2,int xgi,real x,
1620                           int *dx0,int *dx1,real *dcx2)
1621 {
1622     real dcx,tmp;
1623     int  xgi0,xgi1,i;
1624
1625     if (xgi < 0)
1626     {
1627         *dx0 = 0;
1628         xgi0 = -1;
1629         *dx1 = -1;
1630         xgi1 = 0;
1631     }
1632     else if (xgi >= Nx)
1633     {
1634         *dx0 = Nx;
1635         xgi0 = Nx-1;
1636         *dx1 = Nx-1;
1637         xgi1 = Nx;
1638     }
1639     else
1640     {
1641         dcx2[xgi] = 0;
1642         *dx0 = xgi;
1643         xgi0 = xgi-1;
1644         *dx1 = xgi;
1645         xgi1 = xgi+1;
1646     }
1647
1648     for(i=xgi0; i>=0; i--)
1649     {
1650         dcx = (i+1)*gridx-x;
1651         tmp = dcx*dcx;
1652         if (tmp >= rc2)
1653             break;
1654         *dx0 = i;
1655         dcx2[i] = tmp;
1656     }
1657     for(i=xgi1; i<Nx; i++)
1658     {
1659         dcx = i*gridx-x;
1660         tmp = dcx*dcx;
1661         if (tmp >= rc2)
1662         {
1663             break;
1664         }
1665         *dx1 = i;
1666         dcx2[i] = tmp;
1667     }
1668 }
1669
1670 static inline void get_dx_dd(int Nx,real gridx,real rc2,int xgi,real x,
1671                              int ncpddc,int shift_min,int shift_max,
1672                              int *g0,int *g1,real *dcx2)
1673 {
1674     real dcx,tmp;
1675     int  g_min,g_max,shift_home;
1676
1677     if (xgi < 0)
1678     {
1679         g_min = 0;
1680         g_max = Nx - 1;
1681         *g0   = 0;
1682         *g1   = -1;
1683     }
1684     else if (xgi >= Nx)
1685     {
1686         g_min = 0;
1687         g_max = Nx - 1;
1688         *g0   = Nx;
1689         *g1   = Nx - 1;
1690     }
1691     else
1692     {
1693         if (ncpddc == 0)
1694         {
1695             g_min = 0;
1696             g_max = Nx - 1;
1697         }
1698         else
1699         {
1700             if (xgi < ncpddc)
1701             {
1702                 shift_home = 0;
1703             }
1704             else
1705             {
1706                 shift_home = -1;
1707             }
1708             g_min = (shift_min == shift_home ? 0          : ncpddc);
1709             g_max = (shift_max == shift_home ? ncpddc - 1 : Nx - 1);
1710         }
1711         if (shift_min > 0)
1712         {
1713             *g0 = g_min;
1714             *g1 = g_min - 1;
1715         }
1716         else if (shift_max < 0)
1717         {
1718             *g0 = g_max + 1;
1719             *g1 = g_max;
1720         }
1721         else
1722         {
1723             *g0 = xgi;
1724             *g1 = xgi;
1725             dcx2[xgi] = 0;
1726         }
1727     }
1728
1729     while (*g0 > g_min)
1730     {
1731         /* Check one grid cell down */
1732         dcx = ((*g0 - 1) + 1)*gridx - x;
1733         tmp = dcx*dcx;
1734         if (tmp >= rc2)
1735         {
1736             break;
1737         }
1738         (*g0)--;
1739         dcx2[*g0] = tmp;
1740     }
1741
1742     while (*g1 < g_max)
1743     {
1744         /* Check one grid cell up */
1745         dcx = (*g1 + 1)*gridx - x;
1746         tmp = dcx*dcx;
1747         if (tmp >= rc2)
1748         {
1749             break;
1750         }
1751         (*g1)++;
1752         dcx2[*g1] = tmp;
1753     }
1754 }
1755
1756
1757 #define sqr(x) ((x)*(x))
1758 #define calc_dx2(XI,YI,ZI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]) + sqr(ZI-y[ZZ]))
1759 #define calc_cyl_dx2(XI,YI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]))
1760 /****************************************************
1761  *
1762  *    F A S T   N E I G H B O R  S E A R C H I N G
1763  *
1764  *    Optimized neighboursearching routine using grid
1765  *    at least 1x1x1, see GROMACS manual
1766  *
1767  ****************************************************/
1768
1769 static void do_longrange(t_commrec *cr,gmx_localtop_t *top,t_forcerec *fr,
1770                          int ngid,t_mdatoms *md,int icg,
1771                          int jgid,int nlr,
1772                          atom_id lr[],t_excl bexcl[],int shift,
1773                          rvec x[],rvec box_size,t_nrnb *nrnb,
1774                          real lambda,real *dvdlambda,
1775                          gmx_grppairener_t *grppener,
1776                          bool bDoVdW,bool bDoCoul,
1777                          bool bEvaluateNow,put_in_list_t *put_in_list,
1778                          bool bHaveVdW[],
1779                          bool bDoForces,rvec *f)
1780 {
1781     int n,i;
1782     t_nblist *nl;
1783
1784     for(n=0; n<fr->nnblists; n++)
1785     {
1786         for(i=0; (i<eNL_NR); i++)
1787         {
1788             nl = &fr->nblists[n].nlist_lr[i];
1789             if ((nl->nri > nl->maxnri-32) || bEvaluateNow)
1790             {
1791                 close_neighbor_list(fr,TRUE,n,i,FALSE);
1792                 /* Evaluate the energies and forces */
1793                 do_nonbonded(cr,fr,x,f,md,NULL,
1794                              grppener->ener[fr->bBHAM ? egBHAMLR : egLJLR],
1795                              grppener->ener[egCOULLR],
1796                                                          grppener->ener[egGB],box_size,
1797                              nrnb,lambda,dvdlambda,n,i,
1798                              GMX_DONB_LR | GMX_DONB_FORCES);
1799
1800                 reset_neighbor_list(fr,TRUE,n,i);
1801             }
1802         }
1803     }
1804
1805     if (!bEvaluateNow)
1806     {
1807         /* Put the long range particles in a list */
1808         /* do_longrange is never called for QMMM  */
1809         put_in_list(bHaveVdW,ngid,md,icg,jgid,nlr,lr,top->cgs.index,
1810                     bexcl,shift,fr,TRUE,bDoVdW,bDoCoul);
1811     }
1812 }
1813
1814 static void get_cutoff2(t_forcerec *fr,bool bDoLongRange,
1815                         real *rvdw2,real *rcoul2,
1816                         real *rs2,real *rm2,real *rl2)
1817 {
1818     *rs2 = sqr(fr->rlist);
1819     if (bDoLongRange && fr->bTwinRange)
1820     {
1821         /* The VdW and elec. LR cut-off's could be different,
1822          * so we can not simply set them to rlistlong.
1823          */
1824         if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(fr->vdwtype) &&
1825             fr->rvdw > fr->rlist)
1826         {
1827             *rvdw2  = sqr(fr->rlistlong);
1828         }
1829         else
1830         {
1831             *rvdw2  = sqr(fr->rvdw);
1832         }
1833         if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(fr->eeltype) &&
1834             fr->rcoulomb > fr->rlist)
1835         {
1836             *rcoul2 = sqr(fr->rlistlong);
1837         }
1838         else
1839         {
1840             *rcoul2 = sqr(fr->rcoulomb);
1841         }
1842     }
1843     else
1844     {
1845         /* Workaround for a gcc -O3 or -ffast-math problem */
1846         *rvdw2  = *rs2;
1847         *rcoul2 = *rs2;
1848     }
1849     *rm2 = min(*rvdw2,*rcoul2);
1850     *rl2 = max(*rvdw2,*rcoul2);
1851 }
1852
1853 static void init_nsgrid_lists(t_forcerec *fr,int ngid,gmx_ns_t *ns)
1854 {
1855     real rvdw2,rcoul2,rs2,rm2,rl2;
1856     int j;
1857
1858     get_cutoff2(fr,TRUE,&rvdw2,&rcoul2,&rs2,&rm2,&rl2);
1859
1860     /* Short range buffers */
1861     snew(ns->nl_sr,ngid);
1862     /* Counters */
1863     snew(ns->nsr,ngid);
1864     snew(ns->nlr_ljc,ngid);
1865     snew(ns->nlr_one,ngid);
1866
1867     if (rm2 > rs2)
1868     {
1869             /* Long range VdW and Coul buffers */
1870         snew(ns->nl_lr_ljc,ngid);
1871     }
1872     if (rl2 > rm2)
1873     {
1874         /* Long range VdW or Coul only buffers */
1875         snew(ns->nl_lr_one,ngid);
1876     }
1877     for(j=0; (j<ngid); j++) {
1878         snew(ns->nl_sr[j],MAX_CG);
1879         if (rm2 > rs2)
1880         {
1881             snew(ns->nl_lr_ljc[j],MAX_CG);
1882         }
1883         if (rl2 > rm2)
1884         {
1885             snew(ns->nl_lr_one[j],MAX_CG);
1886         }
1887     }
1888     if (debug)
1889     {
1890         fprintf(debug,
1891                 "ns5_core: rs2 = %g, rm2 = %g, rl2 = %g (nm^2)\n",
1892                 rs2,rm2,rl2);
1893     }
1894 }
1895
1896 static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
1897                        matrix box,rvec box_size,int ngid,
1898                        gmx_localtop_t *top,
1899                        t_grid *grid,rvec x[],
1900                        t_excl bexcl[],bool *bExcludeAlleg,
1901                        t_nrnb *nrnb,t_mdatoms *md,
1902                        real lambda,real *dvdlambda,
1903                        gmx_grppairener_t *grppener,
1904                        put_in_list_t *put_in_list,
1905                        bool bHaveVdW[],
1906                        bool bDoLongRange,bool bDoForces,rvec *f,
1907                        bool bMakeQMMMnblist)
1908 {
1909     gmx_ns_t *ns;
1910     atom_id **nl_lr_ljc,**nl_lr_one,**nl_sr;
1911     int     *nlr_ljc,*nlr_one,*nsr;
1912     gmx_domdec_t *dd=NULL;
1913     t_block *cgs=&(top->cgs);
1914     int     *cginfo=fr->cginfo;
1915     /* atom_id *i_atoms,*cgsindex=cgs->index; */
1916     ivec    sh0,sh1,shp;
1917     int     cell_x,cell_y,cell_z;
1918     int     d,tx,ty,tz,dx,dy,dz,cj;
1919 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1920     int     zsh_ty,zsh_tx,ysh_tx;
1921 #endif
1922     int     dx0,dx1,dy0,dy1,dz0,dz1;
1923     int     Nx,Ny,Nz,shift=-1,j,nrj,nns,nn=-1;
1924     real    gridx,gridy,gridz,grid_x,grid_y,grid_z;
1925     real    *dcx2,*dcy2,*dcz2;
1926     int     zgi,ygi,xgi;
1927     int     cg0,cg1,icg=-1,cgsnr,i0,igid,nri,naaj,max_jcg;
1928     int     jcg0,jcg1,jjcg,cgj0,jgid;
1929     int     *grida,*gridnra,*gridind;
1930     bool    rvdw_lt_rcoul,rcoul_lt_rvdw;
1931     rvec    xi,*cgcm,grid_offset;
1932     real    r2,rs2,rvdw2,rcoul2,rm2,rl2,XI,YI,ZI,dcx,dcy,dcz,tmp1,tmp2;
1933     int     *i_egp_flags;
1934     bool    bDomDec,bTriclinicX,bTriclinicY;
1935     ivec    ncpddc;
1936
1937     ns = &fr->ns;
1938
1939     bDomDec = DOMAINDECOMP(cr);
1940     if (bDomDec)
1941     {
1942         dd = cr->dd;
1943     }
1944
1945     bTriclinicX = ((YY < grid->npbcdim &&
1946                     (!bDomDec || dd->nc[YY]==1) && box[YY][XX] != 0) ||
1947                    (ZZ < grid->npbcdim &&
1948                     (!bDomDec || dd->nc[ZZ]==1) && box[ZZ][XX] != 0));
1949     bTriclinicY =  (ZZ < grid->npbcdim &&
1950                     (!bDomDec || dd->nc[ZZ]==1) && box[ZZ][YY] != 0);
1951
1952     cgsnr    = cgs->nr;
1953
1954     get_cutoff2(fr,bDoLongRange,&rvdw2,&rcoul2,&rs2,&rm2,&rl2);
1955
1956     rvdw_lt_rcoul = (rvdw2 >= rcoul2);
1957     rcoul_lt_rvdw = (rcoul2 >= rvdw2);
1958
1959     if (bMakeQMMMnblist)
1960     {
1961         rm2 = rl2;
1962         rs2 = rl2;
1963     }
1964
1965     nl_sr     = ns->nl_sr;
1966     nsr       = ns->nsr;
1967     nl_lr_ljc = ns->nl_lr_ljc;
1968     nl_lr_one = ns->nl_lr_one;
1969     nlr_ljc   = ns->nlr_ljc;
1970     nlr_one   = ns->nlr_one;
1971
1972     /* Unpack arrays */
1973     cgcm    = fr->cg_cm;
1974     Nx      = grid->n[XX];
1975     Ny      = grid->n[YY];
1976     Nz      = grid->n[ZZ];
1977     grida   = grid->a;
1978     gridind = grid->index;
1979     gridnra = grid->nra;
1980     nns     = 0;
1981
1982     gridx      = grid->cell_size[XX];
1983     gridy      = grid->cell_size[YY];
1984     gridz      = grid->cell_size[ZZ];
1985     grid_x     = 1/gridx;
1986     grid_y     = 1/gridy;
1987     grid_z     = 1/gridz;
1988     copy_rvec(grid->cell_offset,grid_offset);
1989     copy_ivec(grid->ncpddc,ncpddc);
1990     dcx2       = grid->dcx2;
1991     dcy2       = grid->dcy2;
1992     dcz2       = grid->dcz2;
1993
1994 #ifdef ALLOW_OFFDIAG_LT_HALFDIAG
1995     zsh_ty = floor(-box[ZZ][YY]/box[YY][YY]+0.5);
1996     zsh_tx = floor(-box[ZZ][XX]/box[XX][XX]+0.5);
1997     ysh_tx = floor(-box[YY][XX]/box[XX][XX]+0.5);
1998     if (zsh_tx!=0 && ysh_tx!=0)
1999     {
2000         /* This could happen due to rounding, when both ratios are 0.5 */
2001         ysh_tx = 0;
2002     }
2003 #endif
2004
2005     debug_gmx();
2006
2007     if (fr->n_tpi)
2008     {
2009         /* We only want a list for the test particle */
2010         cg0 = cgsnr - 1;
2011     }
2012     else
2013     {
2014         cg0 = grid->icg0;
2015     }
2016     cg1 = grid->icg1;
2017
2018     /* Set the shift range */
2019     for(d=0; d<DIM; d++)
2020     {
2021         sh0[d] = -1;
2022         sh1[d] = 1;
2023         /* Check if we need periodicity shifts.
2024          * Without PBC or with domain decomposition we don't need them.
2025          */
2026         if (d >= ePBC2npbcdim(fr->ePBC) || (bDomDec && dd->nc[d] > 1))
2027         {
2028             shp[d] = 0;
2029         }
2030         else
2031         {
2032             if (d == XX &&
2033                 box[XX][XX] - fabs(box[YY][XX]) - fabs(box[ZZ][XX]) < sqrt(rl2))
2034             {
2035                 shp[d] = 2;
2036             }
2037             else
2038             {
2039                 shp[d] = 1;
2040             }
2041         }
2042     }
2043
2044     /* Loop over charge groups */
2045     for(icg=cg0; (icg < cg1); icg++)
2046     {
2047         igid = GET_CGINFO_GID(cginfo[icg]);
2048         /* Skip this charge group if all energy groups are excluded! */
2049         if (bExcludeAlleg[igid])
2050         {
2051             continue;
2052         }
2053
2054         i0   = cgs->index[icg];
2055
2056         if (bMakeQMMMnblist)
2057         {
2058             /* Skip this charge group if it is not a QM atom while making a
2059              * QM/MM neighbourlist
2060              */
2061             if (md->bQM[i0]==FALSE)
2062             {
2063                 continue; /* MM particle, go to next particle */
2064             }
2065
2066             /* Compute the number of charge groups that fall within the control
2067              * of this one (icg)
2068              */
2069             naaj    = calc_naaj(icg,cgsnr);
2070             jcg0    = icg;
2071             jcg1    = icg + naaj;
2072             max_jcg = cgsnr;
2073         }
2074         else
2075         {
2076             /* make a normal neighbourlist */
2077
2078             if (bDomDec)
2079             {
2080                 /* Get the j charge-group and dd cell shift ranges */
2081                 dd_get_ns_ranges(cr->dd,icg,&jcg0,&jcg1,sh0,sh1);
2082                 max_jcg = 0;
2083             }
2084             else
2085             {
2086                 /* Compute the number of charge groups that fall within the control
2087                  * of this one (icg)
2088                  */
2089                 naaj = calc_naaj(icg,cgsnr);
2090                 jcg0 = icg;
2091                 jcg1 = icg + naaj;
2092
2093                 if (fr->n_tpi)
2094                 {
2095                     /* The i-particle is awlways the test particle,
2096                      * so we want all j-particles
2097                      */
2098                     max_jcg = cgsnr - 1;
2099                 }
2100                 else
2101                 {
2102                     max_jcg  = jcg1 - cgsnr;
2103                 }
2104             }
2105         }
2106
2107         i_egp_flags = fr->egp_flags + igid*ngid;
2108
2109         /* Set the exclusions for the atoms in charge group icg using a bitmask */
2110         setexcl(i0,cgs->index[icg+1],&top->excls,TRUE,bexcl);
2111
2112         ci2xyz(grid,icg,&cell_x,&cell_y,&cell_z);
2113
2114         /* Changed iicg to icg, DvdS 990115
2115          * (but see consistency check above, DvdS 990330)
2116          */
2117 #ifdef NS5DB
2118         fprintf(log,"icg=%5d, naaj=%5d, cell %d %d %d\n",
2119                 icg,naaj,cell_x,cell_y,cell_z);
2120 #endif
2121         /* Loop over shift vectors in three dimensions */
2122         for (tz=-shp[ZZ]; tz<=shp[ZZ]; tz++)
2123         {
2124             ZI = cgcm[icg][ZZ]+tz*box[ZZ][ZZ];
2125             /* Calculate range of cells in Z direction that have the shift tz */
2126             zgi = cell_z + tz*Nz;
2127 #define FAST_DD_NS
2128 #ifndef FAST_DD_NS
2129             get_dx(Nz,gridz,rl2,zgi,ZI,&dz0,&dz1,dcz2);
2130 #else
2131             get_dx_dd(Nz,gridz,rl2,zgi,ZI-grid_offset[ZZ],
2132                       ncpddc[ZZ],sh0[ZZ],sh1[ZZ],&dz0,&dz1,dcz2);
2133 #endif
2134             if (dz0 > dz1)
2135             {
2136                 continue;
2137             }
2138             for (ty=-shp[YY]; ty<=shp[YY]; ty++)
2139             {
2140                 YI = cgcm[icg][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
2141                 /* Calculate range of cells in Y direction that have the shift ty */
2142                 if (bTriclinicY)
2143                 {
2144                     ygi = (int)(Ny + (YI - grid_offset[YY])*grid_y) - Ny;
2145                 }
2146                 else
2147                 {
2148                     ygi = cell_y + ty*Ny;
2149                 }
2150 #ifndef FAST_DD_NS
2151                 get_dx(Ny,gridy,rl2,ygi,YI,&dy0,&dy1,dcy2);
2152 #else
2153                 get_dx_dd(Ny,gridy,rl2,ygi,YI-grid_offset[YY],
2154                           ncpddc[YY],sh0[YY],sh1[YY],&dy0,&dy1,dcy2);
2155 #endif
2156                 if (dy0 > dy1)
2157                 {
2158                     continue;
2159                 }
2160                 for (tx=-shp[XX]; tx<=shp[XX]; tx++)
2161                 {
2162                     XI = cgcm[icg][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
2163                     /* Calculate range of cells in X direction that have the shift tx */
2164                     if (bTriclinicX)
2165                     {
2166                         xgi = (int)(Nx + (XI - grid_offset[XX])*grid_x) - Nx;
2167                     }
2168                     else
2169                     {
2170                         xgi = cell_x + tx*Nx;
2171                     }
2172 #ifndef FAST_DD_NS
2173                     get_dx(Nx,gridx,rl2,xgi*Nx,XI,&dx0,&dx1,dcx2);
2174 #else
2175                     get_dx_dd(Nx,gridx,rl2,xgi,XI-grid_offset[XX],
2176                               ncpddc[XX],sh0[XX],sh1[XX],&dx0,&dx1,dcx2);
2177 #endif
2178                     if (dx0 > dx1)
2179                     {
2180                         continue;
2181                     }
2182                     /* Get shift vector */
2183                     shift=XYZ2IS(tx,ty,tz);
2184 #ifdef NS5DB
2185                     range_check(shift,0,SHIFTS);
2186 #endif
2187                     for(nn=0; (nn<ngid); nn++)
2188                     {
2189                         nsr[nn]      = 0;
2190                         nlr_ljc[nn]  = 0;
2191                         nlr_one[nn] = 0;
2192                     }
2193 #ifdef NS5DB
2194                     fprintf(log,"shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
2195                             shift,dx0,dx1,dy0,dy1,dz0,dz1);
2196                     fprintf(log,"cgcm: %8.3f  %8.3f  %8.3f\n",cgcm[icg][XX],
2197                             cgcm[icg][YY],cgcm[icg][ZZ]);
2198                     fprintf(log,"xi:   %8.3f  %8.3f  %8.3f\n",XI,YI,ZI);
2199 #endif
2200                     for (dx=dx0; (dx<=dx1); dx++)
2201                     {
2202                         tmp1 = rl2 - dcx2[dx];
2203                         for (dy=dy0; (dy<=dy1); dy++)
2204                         {
2205                             tmp2 = tmp1 - dcy2[dy];
2206                             if (tmp2 > 0)
2207                             {
2208                                 for (dz=dz0; (dz<=dz1); dz++) {
2209                                     if (tmp2 > dcz2[dz]) {
2210                                         /* Find grid-cell cj in which possible neighbours are */
2211                                         cj   = xyz2ci(Ny,Nz,dx,dy,dz);
2212
2213                                         /* Check out how many cgs (nrj) there in this cell */
2214                                         nrj  = gridnra[cj];
2215
2216                                         /* Find the offset in the cg list */
2217                                         cgj0 = gridind[cj];
2218
2219                                         /* Check if all j's are out of range so we
2220                                          * can skip the whole cell.
2221                                          * Should save some time, especially with DD.
2222                                          */
2223                                         if (nrj == 0 ||
2224                                             (grida[cgj0] >= max_jcg &&
2225                                              (grida[cgj0] >= jcg1 || grida[cgj0+nrj-1] < jcg0)))
2226                                         {
2227                                             continue;
2228                                         }
2229
2230                                         /* Loop over cgs */
2231                                         for (j=0; (j<nrj); j++)
2232                                         {
2233                                             jjcg = grida[cgj0+j];
2234
2235                                             /* check whether this guy is in range! */
2236                                             if ((jjcg >= jcg0 && jjcg < jcg1) ||
2237                                                 (jjcg < max_jcg))
2238                                             {
2239                                                 r2=calc_dx2(XI,YI,ZI,cgcm[jjcg]);
2240                                                 if (r2 < rl2) {
2241                                                     /* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
2242                                                     jgid = GET_CGINFO_GID(cginfo[jjcg]);
2243                                                     /* check energy group exclusions */
2244                                                     if (!(i_egp_flags[jgid] & EGP_EXCL))
2245                                                     {
2246                                                         if (r2 < rs2)
2247                                                         {
2248                                                             if (nsr[jgid] >= MAX_CG)
2249                                                             {
2250                                                                 put_in_list(bHaveVdW,ngid,md,icg,jgid,
2251                                                                             nsr[jgid],nl_sr[jgid],
2252                                                                             cgs->index,/* cgsatoms, */ bexcl,
2253                                                                             shift,fr,FALSE,TRUE,TRUE);
2254                                                                 nsr[jgid]=0;
2255                                                             }
2256                                                             nl_sr[jgid][nsr[jgid]++]=jjcg;
2257                                                         }
2258                                                         else if (r2 < rm2)
2259                                                         {
2260                                                             if (nlr_ljc[jgid] >= MAX_CG)
2261                                                             {
2262                                                                 do_longrange(cr,top,fr,ngid,md,icg,jgid,
2263                                                                              nlr_ljc[jgid],
2264                                                                              nl_lr_ljc[jgid],bexcl,shift,x,
2265                                                                              box_size,nrnb,
2266                                                                              lambda,dvdlambda,
2267                                                                              grppener,
2268                                                                              TRUE,TRUE,FALSE,
2269                                                                              put_in_list,
2270                                                                              bHaveVdW,
2271                                                                              bDoForces,f);
2272                                                                 nlr_ljc[jgid]=0;
2273                                                             }
2274                                                             nl_lr_ljc[jgid][nlr_ljc[jgid]++]=jjcg;
2275                                                         }
2276                                                         else
2277                                                         {
2278                                                             if (nlr_one[jgid] >= MAX_CG) {
2279                                                                 do_longrange(cr,top,fr,ngid,md,icg,jgid,
2280                                                                              nlr_one[jgid],
2281                                                                              nl_lr_one[jgid],bexcl,shift,x,
2282                                                                              box_size,nrnb,
2283                                                                              lambda,dvdlambda,
2284                                                                              grppener,
2285                                                                              rvdw_lt_rcoul,rcoul_lt_rvdw,FALSE,
2286                                                                              put_in_list,
2287                                                                              bHaveVdW,
2288                                                                              bDoForces,f);
2289                                                                 nlr_one[jgid]=0;
2290                                                             }
2291                                                             nl_lr_one[jgid][nlr_one[jgid]++]=jjcg;
2292                                                         }
2293                                                     }
2294                                                 }
2295                                                 nns++;
2296                                             }
2297                                         }
2298                                     }
2299                                 }
2300                             }
2301                         }
2302                     }
2303                     /* CHECK whether there is anything left in the buffers */
2304                     for(nn=0; (nn<ngid); nn++)
2305                     {
2306                         if (nsr[nn] > 0)
2307                         {
2308                             put_in_list(bHaveVdW,ngid,md,icg,nn,nsr[nn],nl_sr[nn],
2309                                         cgs->index, /* cgsatoms, */ bexcl,
2310                                         shift,fr,FALSE,TRUE,TRUE);
2311                         }
2312
2313                         if (nlr_ljc[nn] > 0)
2314                         {
2315                             do_longrange(cr,top,fr,ngid,md,icg,nn,nlr_ljc[nn],
2316                                          nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
2317                                          lambda,dvdlambda,grppener,TRUE,TRUE,FALSE,
2318                                          put_in_list,bHaveVdW,bDoForces,f);
2319                         }
2320
2321                         if (nlr_one[nn] > 0)
2322                         {
2323                             do_longrange(cr,top,fr,ngid,md,icg,nn,nlr_one[nn],
2324                                          nl_lr_one[nn],bexcl,shift,x,box_size,nrnb,
2325                                          lambda,dvdlambda,grppener,
2326                                          rvdw_lt_rcoul,rcoul_lt_rvdw,FALSE,
2327                                          put_in_list,bHaveVdW,bDoForces,f);
2328                         }
2329                     }
2330                 }
2331             }
2332         }
2333         /* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */
2334         setexcl(cgs->index[icg],cgs->index[icg+1],&top->excls,FALSE,bexcl);
2335     }
2336     /* Perform any left over force calculations */
2337     for (nn=0; (nn<ngid); nn++)
2338     {
2339         if (rm2 > rs2)
2340         {
2341             do_longrange(cr,top,fr,0,md,icg,nn,nlr_ljc[nn],
2342                          nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
2343                          lambda,dvdlambda,grppener,
2344                          TRUE,TRUE,TRUE,put_in_list,bHaveVdW,bDoForces,f);
2345         }
2346         if (rl2 > rm2) {
2347             do_longrange(cr,top,fr,0,md,icg,nn,nlr_one[nn],
2348                          nl_lr_one[nn],bexcl,shift,x,box_size,nrnb,
2349                          lambda,dvdlambda,grppener,
2350                          rvdw_lt_rcoul,rcoul_lt_rvdw,
2351                          TRUE,put_in_list,bHaveVdW,bDoForces,f);
2352         }
2353     }
2354     debug_gmx();
2355
2356     /* Close off short range neighbourlists */
2357     close_neighbor_list(fr,FALSE,-1,-1,bMakeQMMMnblist);
2358
2359     return nns;
2360 }
2361
2362 void ns_realloc_natoms(gmx_ns_t *ns,int natoms)
2363 {
2364     int i;
2365
2366     if (natoms > ns->nra_alloc)
2367     {
2368         ns->nra_alloc = over_alloc_dd(natoms);
2369         srenew(ns->bexcl,ns->nra_alloc);
2370         for(i=0; i<ns->nra_alloc; i++)
2371         {
2372             ns->bexcl[i] = 0;
2373         }
2374     }
2375 }
2376
2377 void init_ns(FILE *fplog,const t_commrec *cr,
2378              gmx_ns_t *ns,t_forcerec *fr,
2379              const gmx_mtop_t *mtop,
2380              matrix box)
2381 {
2382     int  mt,icg,nr_in_cg,maxcg,i,j,jcg,ngid,ncg;
2383     t_block *cgs;
2384     char *ptr;
2385
2386     /* Compute largest charge groups size (# atoms) */
2387     nr_in_cg=1;
2388     for(mt=0; mt<mtop->nmoltype; mt++) {
2389         cgs = &mtop->moltype[mt].cgs;
2390         for (icg=0; (icg < cgs->nr); icg++)
2391         {
2392             nr_in_cg=max(nr_in_cg,(int)(cgs->index[icg+1]-cgs->index[icg]));
2393         }
2394     }
2395
2396     /* Verify whether largest charge group is <= max cg.
2397      * This is determined by the type of the local exclusion type
2398      * Exclusions are stored in bits. (If the type is not large
2399      * enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
2400      */
2401     maxcg = sizeof(t_excl)*8;
2402     if (nr_in_cg > maxcg)
2403     {
2404         gmx_fatal(FARGS,"Max #atoms in a charge group: %d > %d\n",
2405                   nr_in_cg,maxcg);
2406     }
2407
2408     ngid = mtop->groups.grps[egcENER].nr;
2409     snew(ns->bExcludeAlleg,ngid);
2410     for(i=0; i<ngid; i++) {
2411         ns->bExcludeAlleg[i] = TRUE;
2412         for(j=0; j<ngid; j++)
2413         {
2414             if (!(fr->egp_flags[i*ngid+j] & EGP_EXCL))
2415             {
2416                 ns->bExcludeAlleg[i] = FALSE;
2417             }
2418         }
2419     }
2420
2421     if (fr->bGrid) {
2422         /* Grid search */
2423         ns->grid = init_grid(fplog,fr);
2424         init_nsgrid_lists(fr,ngid,ns);
2425     }
2426     else
2427     {
2428         /* Simple search */
2429         snew(ns->ns_buf,ngid);
2430         for(i=0; (i<ngid); i++)
2431         {
2432             snew(ns->ns_buf[i],SHIFTS);
2433         }
2434         ncg = ncg_mtop(mtop);
2435         snew(ns->simple_aaj,2*ncg);
2436         for(jcg=0; (jcg<ncg); jcg++)
2437         {
2438             ns->simple_aaj[jcg]     = jcg;
2439             ns->simple_aaj[jcg+ncg] = jcg;
2440         }
2441     }
2442
2443     /* Create array that determines whether or not atoms have VdW */
2444     snew(ns->bHaveVdW,fr->ntype);
2445     for(i=0; (i<fr->ntype); i++)
2446     {
2447         for(j=0; (j<fr->ntype); j++)
2448         {
2449             ns->bHaveVdW[i] = (ns->bHaveVdW[i] ||
2450                                (fr->bBHAM ?
2451                                 ((BHAMA(fr->nbfp,fr->ntype,i,j) != 0) ||
2452                                  (BHAMB(fr->nbfp,fr->ntype,i,j) != 0) ||
2453                                  (BHAMC(fr->nbfp,fr->ntype,i,j) != 0)) :
2454                                 ((C6(fr->nbfp,fr->ntype,i,j) != 0) ||
2455                                  (C12(fr->nbfp,fr->ntype,i,j) != 0))));
2456         }
2457     }
2458     if (debug)
2459         pr_bvec(debug,0,"bHaveVdW",ns->bHaveVdW,fr->ntype,TRUE);
2460
2461     ns->nra_alloc = 0;
2462     ns->bexcl = NULL;
2463     if (!DOMAINDECOMP(cr))
2464     {
2465         /* This could be reduced with particle decomposition */
2466         ns_realloc_natoms(ns,mtop->natoms);
2467     }
2468
2469     ns->nblist_initialized=FALSE;
2470
2471     /* nbr list debug dump */
2472     {
2473         char *ptr=getenv("GMX_DUMP_NL");
2474         if (ptr)
2475         {
2476             ns->dump_nl=strtol(ptr,NULL,10);
2477             if (fplog)
2478             {
2479                 fprintf(fplog, "GMX_DUMP_NL = %d", ns->dump_nl);
2480             }
2481         }
2482         else
2483         {
2484             ns->dump_nl=0;
2485         }
2486     }
2487 }
2488
2489
2490 int search_neighbours(FILE *log,t_forcerec *fr,
2491                       rvec x[],matrix box,
2492                       gmx_localtop_t *top,
2493                       gmx_groups_t *groups,
2494                       t_commrec *cr,
2495                       t_nrnb *nrnb,t_mdatoms *md,
2496                       real lambda,real *dvdlambda,
2497                       gmx_grppairener_t *grppener,
2498                       bool bFillGrid,
2499                       bool bDoLongRange,
2500                       bool bDoForces,rvec *f)
2501 {
2502     t_block  *cgs=&(top->cgs);
2503     rvec     box_size,grid_x0,grid_x1;
2504     int      i,j,m,ngid;
2505     real     min_size,grid_dens;
2506     int      nsearch;
2507     bool     bGrid;
2508     char     *ptr;
2509     bool     *i_egp_flags;
2510     int      cg_start,cg_end,start,end;
2511     gmx_ns_t *ns;
2512     t_grid   *grid;
2513     gmx_domdec_zones_t *dd_zones;
2514     put_in_list_t *put_in_list;
2515
2516     ns = &fr->ns;
2517
2518     /* Set some local variables */
2519     bGrid = fr->bGrid;
2520     ngid = groups->grps[egcENER].nr;
2521
2522     for(m=0; (m<DIM); m++)
2523     {
2524         box_size[m] = box[m][m];
2525     }
2526
2527     if (fr->ePBC != epbcNONE)
2528     {
2529         if (sqr(fr->rlistlong) >= max_cutoff2(fr->ePBC,box))
2530         {
2531             gmx_fatal(FARGS,"One of the box vectors has become shorter than twice the cut-off length or box_yy-|box_zy| or box_zz has become smaller than the cut-off.");
2532         }
2533         if (!bGrid)
2534         {
2535             min_size = min(box_size[XX],min(box_size[YY],box_size[ZZ]));
2536             if (2*fr->rlistlong >= min_size)
2537                 gmx_fatal(FARGS,"One of the box diagonal elements has become smaller than twice the cut-off length.");
2538         }
2539     }
2540
2541     if (DOMAINDECOMP(cr))
2542     {
2543         ns_realloc_natoms(ns,cgs->index[cgs->nr]);
2544     }
2545     debug_gmx();
2546
2547     /* Reset the neighbourlists */
2548     reset_neighbor_list(fr,FALSE,-1,-1);
2549
2550     if (bGrid && bFillGrid)
2551     {
2552
2553         grid = ns->grid;
2554         if (DOMAINDECOMP(cr))
2555         {
2556             dd_zones = domdec_zones(cr->dd);
2557         }
2558         else
2559         {
2560             dd_zones = NULL;
2561
2562             get_nsgrid_boundaries(grid,NULL,box,NULL,NULL,NULL,
2563                                   cgs->nr,fr->cg_cm,grid_x0,grid_x1,&grid_dens);
2564
2565             grid_first(log,grid,NULL,NULL,fr->ePBC,box,grid_x0,grid_x1,
2566                        fr->rlistlong,grid_dens);
2567         }
2568         debug_gmx();
2569
2570         /* Don't know why this all is... (DvdS 3/99) */
2571 #ifndef SEGV
2572         start = 0;
2573         end   = cgs->nr;
2574 #else
2575         start = fr->cg0;
2576         end   = (cgs->nr+1)/2;
2577 #endif
2578
2579         if (DOMAINDECOMP(cr))
2580         {
2581             end = cgs->nr;
2582             fill_grid(log,dd_zones,grid,end,-1,end,fr->cg_cm);
2583             grid->icg0 = 0;
2584             grid->icg1 = dd_zones->izone[dd_zones->nizone-1].cg1;
2585         }
2586         else
2587         {
2588             fill_grid(log,NULL,grid,cgs->nr,fr->cg0,fr->hcg,fr->cg_cm);
2589             grid->icg0 = fr->cg0;
2590             grid->icg1 = fr->hcg;
2591             debug_gmx();
2592
2593             if (PARTDECOMP(cr))
2594                 mv_grid(cr,grid);
2595             debug_gmx();
2596         }
2597
2598         calc_elemnr(log,grid,start,end,cgs->nr);
2599         calc_ptrs(grid);
2600         grid_last(log,grid,start,end,cgs->nr);
2601
2602         if (gmx_debug_at)
2603         {
2604             check_grid(debug,grid);
2605             print_grid(debug,grid);
2606         }
2607     }
2608     else if (fr->n_tpi)
2609     {
2610         /* Set the grid cell index for the test particle only.
2611          * The cell to cg index is not corrected, but that does not matter.
2612          */
2613         fill_grid(log,NULL,ns->grid,fr->hcg,fr->hcg-1,fr->hcg,fr->cg_cm);
2614     }
2615     debug_gmx();
2616
2617     if (!fr->ns.bCGlist)
2618     {
2619         put_in_list = put_in_list_at;
2620     }
2621     else
2622     {
2623         put_in_list = put_in_list_cg;
2624     }
2625
2626     /* Do the core! */
2627     if (bGrid)
2628     {
2629         grid = ns->grid;
2630         nsearch = nsgrid_core(log,cr,fr,box,box_size,ngid,top,
2631                               grid,x,ns->bexcl,ns->bExcludeAlleg,
2632                               nrnb,md,lambda,dvdlambda,grppener,
2633                               put_in_list,ns->bHaveVdW,
2634                               bDoLongRange,bDoForces,f,
2635                               FALSE);
2636
2637         /* neighbour searching withouth QMMM! QM atoms have zero charge in
2638          * the classical calculation. The charge-charge interaction
2639          * between QM and MM atoms is handled in the QMMM core calculation
2640          * (see QMMM.c). The VDW however, we'd like to compute classically
2641          * and the QM MM atom pairs have just been put in the
2642          * corresponding neighbourlists. in case of QMMM we still need to
2643          * fill a special QMMM neighbourlist that contains all neighbours
2644          * of the QM atoms. If bQMMM is true, this list will now be made:
2645          */
2646         if (fr->bQMMM && fr->qr->QMMMscheme!=eQMMMschemeoniom)
2647         {
2648             nsearch += nsgrid_core(log,cr,fr,box,box_size,ngid,top,
2649                                    grid,x,ns->bexcl,ns->bExcludeAlleg,
2650                                    nrnb,md,lambda,dvdlambda,grppener,
2651                                    put_in_list_qmmm,ns->bHaveVdW,
2652                                    bDoLongRange,bDoForces,f,
2653                                    TRUE);
2654         }
2655     }
2656     else
2657     {
2658         nsearch = ns_simple_core(fr,top,md,box,box_size,
2659                                  ns->bexcl,ns->simple_aaj,
2660                                  ngid,ns->ns_buf,put_in_list,ns->bHaveVdW);
2661     }
2662     debug_gmx();
2663
2664 #ifdef DEBUG
2665     pr_nsblock(log);
2666 #endif
2667
2668     inc_nrnb(nrnb,eNR_NS,nsearch);
2669     /* inc_nrnb(nrnb,eNR_LR,fr->nlr); */
2670
2671     return nsearch;
2672 }
2673
2674 int natoms_beyond_ns_buffer(t_inputrec *ir,t_forcerec *fr,t_block *cgs,
2675                             matrix scale_tot,rvec *x)
2676 {
2677     int  cg0,cg1,cg,a0,a1,a,i,j;
2678     real rint,hbuf2,scale;
2679     rvec *cg_cm,cgsc;
2680     bool bIsotropic;
2681     int  nBeyond;
2682
2683     nBeyond = 0;
2684
2685     rint = max(ir->rcoulomb,ir->rvdw);
2686     if (ir->rlist < rint)
2687     {
2688         gmx_fatal(FARGS,"The neighbor search buffer has negative size: %f nm",
2689                   ir->rlist - rint);
2690     }
2691     cg_cm = fr->cg_cm;
2692
2693     cg0 = fr->cg0;
2694     cg1 = fr->hcg;
2695
2696     if (!EI_DYNAMICS(ir->eI) || !DYNAMIC_BOX(*ir))
2697     {
2698         hbuf2 = sqr(0.5*(ir->rlist - rint));
2699         for(cg=cg0; cg<cg1; cg++)
2700         {
2701             a0 = cgs->index[cg];
2702             a1 = cgs->index[cg+1];
2703             for(a=a0; a<a1; a++)
2704             {
2705                 if (distance2(cg_cm[cg],x[a]) > hbuf2)
2706                 {
2707                     nBeyond++;
2708                 }
2709             }
2710         }
2711     }
2712     else
2713     {
2714         bIsotropic = TRUE;
2715         scale = scale_tot[0][0];
2716         for(i=1; i<DIM; i++)
2717         {
2718             /* With anisotropic scaling, the original spherical ns volumes become
2719              * ellipsoids. To avoid costly transformations we use the minimum
2720              * eigenvalue of the scaling matrix for determining the buffer size.
2721              * Since the lower half is 0, the eigenvalues are the diagonal elements.
2722              */
2723             scale = min(scale,scale_tot[i][i]);
2724             if (scale_tot[i][i] != scale_tot[i-1][i-1])
2725             {
2726                 bIsotropic = FALSE;
2727             }
2728             for(j=0; j<i; j++)
2729             {
2730                 if (scale_tot[i][j] != 0)
2731                 {
2732                     bIsotropic = FALSE;
2733                 }
2734             }
2735         }
2736         hbuf2 = sqr(0.5*(scale*ir->rlist - rint));
2737         if (bIsotropic)
2738         {
2739             for(cg=cg0; cg<cg1; cg++)
2740             {
2741                 svmul(scale,cg_cm[cg],cgsc);
2742                 a0 = cgs->index[cg];
2743                 a1 = cgs->index[cg+1];
2744                 for(a=a0; a<a1; a++)
2745                 {
2746                     if (distance2(cgsc,x[a]) > hbuf2)
2747                     {
2748                         nBeyond++;
2749                     }
2750                 }
2751             }
2752         }
2753         else
2754         {
2755             /* Anistropic scaling */
2756             for(cg=cg0; cg<cg1; cg++)
2757             {
2758                 /* Since scale_tot contains the transpose of the scaling matrix,
2759                  * we need to multiply with the transpose.
2760                  */
2761                 tmvmul_ur0(scale_tot,cg_cm[cg],cgsc);
2762                 a0 = cgs->index[cg];
2763                 a1 = cgs->index[cg+1];
2764                 for(a=a0; a<a1; a++)
2765                 {
2766                     if (distance2(cgsc,x[a]) > hbuf2)
2767                     {
2768                         nBeyond++;
2769                     }
2770                 }
2771             }
2772         }
2773     }
2774
2775     return nBeyond;
2776 }