src/gromacs/mdlib/ns.c

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