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