src/gromacs/mdlib/domdec_top.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2006,2007,2008,2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 #include "gmxpre.h"
  37
  38 #include <string.h>
  39
  40 #include "gromacs/legacyheaders/chargegroup.h"
  41 #include "gromacs/legacyheaders/domdec.h"
  42 #include "gromacs/legacyheaders/domdec_network.h"
  43 #include "gromacs/legacyheaders/force.h"
  44 #include "gromacs/legacyheaders/gmx_ga2la.h"
  45 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
  46 #include "gromacs/legacyheaders/names.h"
  47 #include "gromacs/legacyheaders/network.h"
  48 #include "gromacs/legacyheaders/typedefs.h"
  49 #include "gromacs/legacyheaders/vsite.h"
  50 #include "gromacs/legacyheaders/types/commrec.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/pbcutil/mshift.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/topology/mtop_util.h"
  55 #include "gromacs/topology/topsort.h"
  56 #include "gromacs/utility/cstringutil.h"
  57 #include "gromacs/utility/fatalerror.h"
  58 #include "gromacs/utility/gmxassert.h"
  59 #include "gromacs/utility/smalloc.h"
  60
  61 /* for dd_init_local_state */
  62 #define NITEM_DD_INIT_LOCAL_STATE 5
  63
  64 typedef struct {
  65     int  *index;  /* Index for each atom into il                  */
  66     int  *il;     /* ftype|type|a0|...|an|ftype|...               */
  67 } gmx_reverse_ilist_t;
  68
  69 typedef struct {
  70     int  a_start;
  71     int  a_end;
  72     int  natoms_mol;
  73     int  type;
  74 } gmx_molblock_ind_t;
  75
  76 typedef struct gmx_reverse_top {
  77     gmx_bool             bExclRequired; /* Do we require all exclusions to be assigned? */
  78     gmx_bool             bConstr;       /* Are there constraints in this revserse top?  */
  79     gmx_bool             bSettle;       /* Are there settles in this revserse top?  */
  80     gmx_bool             bBCheck;       /* All bonded interactions have to be assigned? */
  81     gmx_bool             bMultiCGmols;  /* Are the multi charge-group molecules?        */
  82     gmx_reverse_ilist_t *ril_mt;        /* Reverse ilist for all moltypes      */
  83     int                  ril_mt_tot_size;
  84     int                  ilsort;        /* The sorting state of bondeds for free energy */
  85     gmx_molblock_ind_t  *mbi;
  86     int                  nmolblock;
  87
  88     /* Work data structures for multi-threading */
  89     int         nthread;
  90     t_idef     *idef_thread;
  91     int      ***vsite_pbc;
  92     int       **vsite_pbc_nalloc;
  93     int        *nbonded_thread;
  94     t_blocka   *excl_thread;
  95     int        *excl_count_thread;
  96
  97     /* Pointers only used for an error message */
  98     gmx_mtop_t     *err_top_global;
  99     gmx_localtop_t *err_top_local;
 100 } gmx_reverse_top_t;
 101
 102 static int nral_rt(int ftype)
 103 {
 104     /* Returns the number of atom entries for il in gmx_reverse_top_t */
 105     int nral;
 106
 107     nral = NRAL(ftype);
 108     if (interaction_function[ftype].flags & IF_VSITE)
 109     {
 110         /* With vsites the reverse topology contains
 111          * two extra entries for PBC.
 112          */
 113         nral += 2;
 114     }
 115
 116     return nral;
 117 }
 118
 119 /* This function tells which interactions need to be assigned exactly once */
 120 static gmx_bool dd_check_ftype(int ftype, gmx_bool bBCheck,
 121                                gmx_bool bConstr, gmx_bool bSettle)
 122 {
 123     return (((interaction_function[ftype].flags & IF_BOND) &&
 124              !(interaction_function[ftype].flags & IF_VSITE) &&
 125              (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))) ||
 126             (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) ||
 127             (bSettle && ftype == F_SETTLE));
 128 }
 129
 130 static void print_error_header(FILE *fplog, char *moltypename, int nprint)
 131 {
 132     fprintf(fplog, "\nMolecule type '%s'\n", moltypename);
 133     fprintf(stderr, "\nMolecule type '%s'\n", moltypename);
 134     fprintf(fplog,
 135             "the first %d missing interactions, except for exclusions:\n",
 136             nprint);
 137     fprintf(stderr,
 138             "the first %d missing interactions, except for exclusions:\n",
 139             nprint);
 140 }
 141
 142 static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
 143                                           gmx_reverse_top_t *rt,
 144                                           char *moltypename,
 145                                           gmx_reverse_ilist_t *ril,
 146                                           int a_start, int a_end,
 147                                           int nat_mol, int nmol,
 148                                           t_idef *idef)
 149 {
 150     int      nril_mol, *assigned, *gatindex;
 151     int      ftype, ftype_j, nral, i, j_mol, j, a0, a0_mol, mol, a;
 152     int      nprint;
 153     t_ilist *il;
 154     t_iatom *ia;
 155     gmx_bool bFound;
 156
 157     nril_mol = ril->index[nat_mol];
 158     snew(assigned, nmol*nril_mol);
 159
 160     gatindex = cr->dd->gatindex;
 161     for (ftype = 0; ftype < F_NRE; ftype++)
 162     {
 163         if (dd_check_ftype(ftype, rt->bBCheck, rt->bConstr, rt->bSettle))
 164         {
 165             nral = NRAL(ftype);
 166             il   = &idef->il[ftype];
 167             ia   = il->iatoms;
 168             for (i = 0; i < il->nr; i += 1+nral)
 169             {
 170                 a0     = gatindex[ia[1]];
 171                 /* Check if this interaction is in
 172                  * the currently checked molblock.
 173                  */
 174                 if (a0 >= a_start && a0 < a_end)
 175                 {
 176                     mol    = (a0 - a_start)/nat_mol;
 177                     a0_mol = (a0 - a_start) - mol*nat_mol;
 178                     j_mol  = ril->index[a0_mol];
 179                     bFound = FALSE;
 180                     while (j_mol < ril->index[a0_mol+1] && !bFound)
 181                     {
 182                         j       = mol*nril_mol + j_mol;
 183                         ftype_j = ril->il[j_mol];
 184                         /* Here we need to check if this interaction has
 185                          * not already been assigned, since we could have
 186                          * multiply defined interactions.
 187                          */
 188                         if (ftype == ftype_j && ia[0] == ril->il[j_mol+1] &&
 189                             assigned[j] == 0)
 190                         {
 191                             /* Check the atoms */
 192                             bFound = TRUE;
 193                             for (a = 0; a < nral; a++)
 194                             {
 195                                 if (gatindex[ia[1+a]] !=
 196                                     a_start + mol*nat_mol + ril->il[j_mol+2+a])
 197                                 {
 198                                     bFound = FALSE;
 199                                 }
 200                             }
 201                             if (bFound)
 202                             {
 203                                 assigned[j] = 1;
 204                             }
 205                         }
 206                         j_mol += 2 + nral_rt(ftype_j);
 207                     }
 208                     if (!bFound)
 209                     {
 210                         gmx_incons("Some interactions seem to be assigned multiple times");
 211                     }
 212                 }
 213                 ia += 1 + nral;
 214             }
 215         }
 216     }
 217
 218     gmx_sumi(nmol*nril_mol, assigned, cr);
 219
 220     nprint = 10;
 221     i      = 0;
 222     for (mol = 0; mol < nmol; mol++)
 223     {
 224         j_mol = 0;
 225         while (j_mol < nril_mol)
 226         {
 227             ftype = ril->il[j_mol];
 228             nral  = NRAL(ftype);
 229             j     = mol*nril_mol + j_mol;
 230             if (assigned[j] == 0 &&
 231                 !(interaction_function[ftype].flags & IF_VSITE))
 232             {
 233                 if (DDMASTER(cr->dd))
 234                 {
 235                     if (i == 0)
 236                     {
 237                         print_error_header(fplog, moltypename, nprint);
 238                     }
 239                     fprintf(fplog, "%20s atoms",
 240                             interaction_function[ftype].longname);
 241                     fprintf(stderr, "%20s atoms",
 242                             interaction_function[ftype].longname);
 243                     for (a = 0; a < nral; a++)
 244                     {
 245                         fprintf(fplog, "%5d", ril->il[j_mol+2+a]+1);
 246                         fprintf(stderr, "%5d", ril->il[j_mol+2+a]+1);
 247                     }
 248                     while (a < 4)
 249                     {
 250                         fprintf(fplog, "     ");
 251                         fprintf(stderr, "     ");
 252                         a++;
 253                     }
 254                     fprintf(fplog, " global");
 255                     fprintf(stderr, " global");
 256                     for (a = 0; a < nral; a++)
 257                     {
 258                         fprintf(fplog, "%6d",
 259                                 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
 260                         fprintf(stderr, "%6d",
 261                                 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
 262                     }
 263                     fprintf(fplog, "\n");
 264                     fprintf(stderr, "\n");
 265                 }
 266                 i++;
 267                 if (i >= nprint)
 268                 {
 269                     break;
 270                 }
 271             }
 272             j_mol += 2 + nral_rt(ftype);
 273         }
 274     }
 275
 276     sfree(assigned);
 277 }
 278
 279 static void print_missing_interactions_atoms(FILE *fplog, t_commrec *cr,
 280                                              gmx_mtop_t *mtop, t_idef *idef)
 281 {
 282     int                mb, a_start, a_end;
 283     gmx_molblock_t    *molb;
 284     gmx_reverse_top_t *rt;
 285
 286     rt = cr->dd->reverse_top;
 287
 288     /* Print the atoms in the missing interactions per molblock */
 289     a_end = 0;
 290     for (mb = 0; mb < mtop->nmolblock; mb++)
 291     {
 292         molb    = &mtop->molblock[mb];
 293         a_start = a_end;
 294         a_end   = a_start + molb->nmol*molb->natoms_mol;
 295
 296         print_missing_interactions_mb(fplog, cr, rt,
 297                                       *(mtop->moltype[molb->type].name),
 298                                       &rt->ril_mt[molb->type],
 299                                       a_start, a_end, molb->natoms_mol,
 300                                       molb->nmol,
 301                                       idef);
 302     }
 303 }
 304
 305 void dd_print_missing_interactions(FILE *fplog, t_commrec *cr, int local_count,  gmx_mtop_t *top_global, t_state *state_local)
 306 {
 307     int             ndiff_tot, cl[F_NRE], n, ndiff, rest_global, rest_local;
 308     int             ftype, nral;
 309     char            buf[STRLEN];
 310     gmx_domdec_t   *dd;
 311     gmx_mtop_t     *err_top_global;
 312     gmx_localtop_t *err_top_local;
 313
 314     dd = cr->dd;
 315
 316     err_top_global = dd->reverse_top->err_top_global;
 317     err_top_local  = dd->reverse_top->err_top_local;
 318
 319     if (fplog)
 320     {
 321         fprintf(fplog, "\nNot all bonded interactions have been properly assigned to the domain decomposition cells\n");
 322         fflush(fplog);
 323     }
 324
 325     ndiff_tot = local_count - dd->nbonded_global;
 326
 327     for (ftype = 0; ftype < F_NRE; ftype++)
 328     {
 329         nral      = NRAL(ftype);
 330         cl[ftype] = err_top_local->idef.il[ftype].nr/(1+nral);
 331     }
 332
 333     gmx_sumi(F_NRE, cl, cr);
 334
 335     if (DDMASTER(dd))
 336     {
 337         if (fplog)
 338         {
 339             fprintf(fplog, "\nA list of missing interactions:\n");
 340         }
 341         fprintf(stderr, "\nA list of missing interactions:\n");
 342         rest_global = dd->nbonded_global;
 343         rest_local  = local_count;
 344         for (ftype = 0; ftype < F_NRE; ftype++)
 345         {
 346             /* In the reverse and local top all constraints are merged
 347              * into F_CONSTR. So in the if statement we skip F_CONSTRNC
 348              * and add these constraints when doing F_CONSTR.
 349              */
 350             if (((interaction_function[ftype].flags & IF_BOND) &&
 351                  (dd->reverse_top->bBCheck
 352                   || !(interaction_function[ftype].flags & IF_LIMZERO)))
 353                 || (dd->reverse_top->bConstr && ftype == F_CONSTR)
 354                 || (dd->reverse_top->bSettle && ftype == F_SETTLE))
 355             {
 356                 n    = gmx_mtop_ftype_count(err_top_global, ftype);
 357                 if (ftype == F_CONSTR)
 358                 {
 359                     n += gmx_mtop_ftype_count(err_top_global, F_CONSTRNC);
 360                 }
 361                 ndiff = cl[ftype] - n;
 362                 if (ndiff != 0)
 363                 {
 364                     sprintf(buf, "%20s of %6d missing %6d",
 365                             interaction_function[ftype].longname, n, -ndiff);
 366                     if (fplog)
 367                     {
 368                         fprintf(fplog, "%s\n", buf);
 369                     }
 370                     fprintf(stderr, "%s\n", buf);
 371                 }
 372                 rest_global -= n;
 373                 rest_local  -= cl[ftype];
 374             }
 375         }
 376
 377         ndiff = rest_local - rest_global;
 378         if (ndiff != 0)
 379         {
 380             sprintf(buf, "%20s of %6d missing %6d", "exclusions",
 381                     rest_global, -ndiff);
 382             if (fplog)
 383             {
 384                 fprintf(fplog, "%s\n", buf);
 385             }
 386             fprintf(stderr, "%s\n", buf);
 387         }
 388     }
 389
 390     print_missing_interactions_atoms(fplog, cr, err_top_global,
 391                                      &err_top_local->idef);
 392     write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr,
 393                  -1, state_local->x, state_local->box);
 394     if (DDMASTER(dd))
 395     {
 396         if (ndiff_tot > 0)
 397         {
 398             gmx_incons("One or more interactions were multiple assigned in the domain decompostion");
 399         }
 400         else
 401         {
 402             gmx_fatal(FARGS, "%d of the %d bonded interactions could not be calculated because some atoms involved moved further apart than the multi-body cut-off distance (%g nm) or the two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds also see option -ddcheck", -ndiff_tot, cr->dd->nbonded_global, dd_cutoff_mbody(cr->dd), dd_cutoff_twobody(cr->dd));
 403         }
 404     }
 405 }
 406
 407 static void global_atomnr_to_moltype_ind(gmx_reverse_top_t *rt, int i_gl,
 408                                          int *mb, int *mt, int *mol, int *i_mol)
 409 {
 410     gmx_molblock_ind_t *mbi   = rt->mbi;
 411     int                 start = 0;
 412     int                 end   =  rt->nmolblock; /* exclusive */
 413     int                 mid;
 414
 415     /* binary search for molblock_ind */
 416     while (TRUE)
 417     {
 418         mid = (start+end)/2;
 419         if (i_gl >= mbi[mid].a_end)
 420         {
 421             start = mid+1;
 422         }
 423         else if (i_gl < mbi[mid].a_start)
 424         {
 425             end = mid;
 426         }
 427         else
 428         {
 429             break;
 430         }
 431     }
 432
 433     *mb  = mid;
 434     mbi += mid;
 435
 436     *mt    = mbi->type;
 437     *mol   = (i_gl - mbi->a_start) / mbi->natoms_mol;
 438     *i_mol = (i_gl - mbi->a_start) - (*mol)*mbi->natoms_mol;
 439 }
 440
 441 static int count_excls(t_block *cgs, t_blocka *excls, int *n_intercg_excl)
 442 {
 443     int n, cg, at0, at1, at, excl, atj;
 444
 445     n               = 0;
 446     *n_intercg_excl = 0;
 447     for (cg = 0; cg < cgs->nr; cg++)
 448     {
 449         at0 = cgs->index[cg];
 450         at1 = cgs->index[cg+1];
 451         for (at = at0; at < at1; at++)
 452         {
 453             for (excl = excls->index[at]; excl < excls->index[at+1]; excl++)
 454             {
 455                 atj = excls->a[excl];
 456                 if (atj > at)
 457                 {
 458                     n++;
 459                     if (atj < at0 || atj >= at1)
 460                     {
 461                         (*n_intercg_excl)++;
 462                     }
 463                 }
 464             }
 465         }
 466     }
 467
 468     return n;
 469 }
 470
 471 static int low_make_reverse_ilist(t_ilist *il_mt, t_atom *atom,
 472                                   int **vsite_pbc,
 473                                   int *count,
 474                                   gmx_bool bConstr, gmx_bool bSettle,
 475                                   gmx_bool bBCheck,
 476                                   int *r_index, int *r_il,
 477                                   gmx_bool bLinkToAllAtoms,
 478                                   gmx_bool bAssign)
 479 {
 480     int      ftype, nral, i, j, nlink, link;
 481     t_ilist *il;
 482     t_iatom *ia;
 483     atom_id  a;
 484     int      nint;
 485     gmx_bool bVSite;
 486
 487     nint = 0;
 488     for (ftype = 0; ftype < F_NRE; ftype++)
 489     {
 490         if ((interaction_function[ftype].flags & (IF_BOND | IF_VSITE)) ||
 491             (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) ||
 492             (bSettle && ftype == F_SETTLE))
 493         {
 494             bVSite = (interaction_function[ftype].flags & IF_VSITE);
 495             nral   = NRAL(ftype);
 496             il     = &il_mt[ftype];
 497             for (i = 0; i < il->nr; i += 1+nral)
 498             {
 499                 ia = il->iatoms + i;
 500                 if (bLinkToAllAtoms)
 501                 {
 502                     if (bVSite)
 503                     {
 504                         /* We don't need the virtual sites for the cg-links */
 505                         nlink = 0;
 506                     }
 507                     else
 508                     {
 509                         nlink = nral;
 510                     }
 511                 }
 512                 else
 513                 {
 514                     /* Couple to the first atom in the interaction */
 515                     nlink = 1;
 516                 }
 517                 for (link = 0; link < nlink; link++)
 518                 {
 519                     a = ia[1+link];
 520                     if (bAssign)
 521                     {
 522                         r_il[r_index[a]+count[a]] =
 523                             (ftype == F_CONSTRNC ? F_CONSTR : ftype);
 524                         r_il[r_index[a]+count[a]+1] = ia[0];
 525                         for (j = 1; j < 1+nral; j++)
 526                         {
 527                             /* Store the molecular atom number */
 528                             r_il[r_index[a]+count[a]+1+j] = ia[j];
 529                         }
 530                     }
 531                     if (interaction_function[ftype].flags & IF_VSITE)
 532                     {
 533                         if (bAssign)
 534                         {
 535                             /* Add an entry to iatoms for storing
 536                              * which of the constructing atoms are
 537                              * vsites again.
 538                              */
 539                             r_il[r_index[a]+count[a]+2+nral] = 0;
 540                             for (j = 2; j < 1+nral; j++)
 541                             {
 542                                 if (atom[ia[j]].ptype == eptVSite)
 543                                 {
 544                                     r_il[r_index[a]+count[a]+2+nral] |= (2<<j);
 545                                 }
 546                             }
 547                             /* Store vsite pbc atom in a second extra entry */
 548                             r_il[r_index[a]+count[a]+2+nral+1] =
 549                                 (vsite_pbc ? vsite_pbc[ftype-F_VSITE2][i/(1+nral)] : -2);
 550                         }
 551                     }
 552                     else
 553                     {
 554                         /* We do not count vsites since they are always
 555                          * uniquely assigned and can be assigned
 556                          * to multiple nodes with recursive vsites.
 557                          */
 558                         if (bBCheck ||
 559                             !(interaction_function[ftype].flags & IF_LIMZERO))
 560                         {
 561                             nint++;
 562                         }
 563                     }
 564                     count[a] += 2 + nral_rt(ftype);
 565                 }
 566             }
 567         }
 568     }
 569
 570     return nint;
 571 }
 572
 573 static int make_reverse_ilist(gmx_moltype_t *molt,
 574                               int **vsite_pbc,
 575                               gmx_bool bConstr, gmx_bool bSettle,
 576                               gmx_bool bBCheck,
 577                               gmx_bool bLinkToAllAtoms,
 578                               gmx_reverse_ilist_t *ril_mt)
 579 {
 580     int nat_mt, *count, i, nint_mt;
 581
 582     /* Count the interactions */
 583     nat_mt = molt->atoms.nr;
 584     snew(count, nat_mt);
 585     low_make_reverse_ilist(molt->ilist, molt->atoms.atom, vsite_pbc,
 586                            count,
 587                            bConstr, bSettle, bBCheck, NULL, NULL,
 588                            bLinkToAllAtoms, FALSE);
 589
 590     snew(ril_mt->index, nat_mt+1);
 591     ril_mt->index[0] = 0;
 592     for (i = 0; i < nat_mt; i++)
 593     {
 594         ril_mt->index[i+1] = ril_mt->index[i] + count[i];
 595         count[i]           = 0;
 596     }
 597     snew(ril_mt->il, ril_mt->index[nat_mt]);
 598
 599     /* Store the interactions */
 600     nint_mt =
 601         low_make_reverse_ilist(molt->ilist, molt->atoms.atom, vsite_pbc,
 602                                count,
 603                                bConstr, bSettle, bBCheck,
 604                                ril_mt->index, ril_mt->il,
 605                                bLinkToAllAtoms, TRUE);
 606
 607     sfree(count);
 608
 609     return nint_mt;
 610 }
 611
 612 static void destroy_reverse_ilist(gmx_reverse_ilist_t *ril)
 613 {
 614     sfree(ril->index);
 615     sfree(ril->il);
 616 }
 617
 618 static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop, gmx_bool bFE,
 619                                            int ***vsite_pbc_molt,
 620                                            gmx_bool bConstr, gmx_bool bSettle,
 621                                            gmx_bool bBCheck, int *nint)
 622 {
 623     int                mt, i, mb;
 624     gmx_reverse_top_t *rt;
 625     int               *nint_mt;
 626     gmx_moltype_t     *molt;
 627     int                thread;
 628
 629     snew(rt, 1);
 630
 631     /* Should we include constraints (for SHAKE) in rt? */
 632     rt->bConstr = bConstr;
 633     rt->bSettle = bSettle;
 634     rt->bBCheck = bBCheck;
 635
 636     rt->bMultiCGmols = FALSE;
 637     snew(nint_mt, mtop->nmoltype);
 638     snew(rt->ril_mt, mtop->nmoltype);
 639     rt->ril_mt_tot_size = 0;
 640     for (mt = 0; mt < mtop->nmoltype; mt++)
 641     {
 642         molt = &mtop->moltype[mt];
 643         if (molt->cgs.nr > 1)
 644         {
 645             rt->bMultiCGmols = TRUE;
 646         }
 647
 648         /* Make the atom to interaction list for this molecule type */
 649         nint_mt[mt] =
 650             make_reverse_ilist(molt, vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
 651                                rt->bConstr, rt->bSettle, rt->bBCheck, FALSE,
 652                                &rt->ril_mt[mt]);
 653
 654         rt->ril_mt_tot_size += rt->ril_mt[mt].index[molt->atoms.nr];
 655     }
 656     if (debug)
 657     {
 658         fprintf(debug, "The total size of the atom to interaction index is %d integers\n", rt->ril_mt_tot_size);
 659     }
 660
 661     *nint = 0;
 662     for (mb = 0; mb < mtop->nmolblock; mb++)
 663     {
 664         *nint += mtop->molblock[mb].nmol*nint_mt[mtop->molblock[mb].type];
 665     }
 666     sfree(nint_mt);
 667
 668     if (bFE && gmx_mtop_bondeds_free_energy(mtop))
 669     {
 670         rt->ilsort = ilsortFE_UNSORTED;
 671     }
 672     else
 673     {
 674         rt->ilsort = ilsortNO_FE;
 675     }
 676
 677     /* Make a molblock index for fast searching */
 678     snew(rt->mbi, mtop->nmolblock);
 679     rt->nmolblock = mtop->nmolblock;
 680     i             = 0;
 681     for (mb = 0; mb < mtop->nmolblock; mb++)
 682     {
 683         rt->mbi[mb].a_start    = i;
 684         i += mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
 685         rt->mbi[mb].a_end      = i;
 686         rt->mbi[mb].natoms_mol = mtop->molblock[mb].natoms_mol;
 687         rt->mbi[mb].type       = mtop->molblock[mb].type;
 688     }
 689
 690     rt->nthread = gmx_omp_nthreads_get(emntDomdec);
 691     snew(rt->idef_thread, rt->nthread);
 692     if (vsite_pbc_molt != NULL)
 693     {
 694         snew(rt->vsite_pbc, rt->nthread);
 695         snew(rt->vsite_pbc_nalloc, rt->nthread);
 696         for (thread = 0; thread < rt->nthread; thread++)
 697         {
 698             snew(rt->vsite_pbc[thread], F_VSITEN-F_VSITE2+1);
 699             snew(rt->vsite_pbc_nalloc[thread], F_VSITEN-F_VSITE2+1);
 700         }
 701     }
 702     snew(rt->nbonded_thread, rt->nthread);
 703     snew(rt->excl_thread, rt->nthread);
 704     snew(rt->excl_count_thread, rt->nthread);
 705
 706     return rt;
 707 }
 708
 709 void dd_make_reverse_top(FILE *fplog,
 710                          gmx_domdec_t *dd, gmx_mtop_t *mtop,
 711                          gmx_vsite_t *vsite,
 712                          t_inputrec *ir, gmx_bool bBCheck)
 713 {
 714     int             mb, nexcl, nexcl_icg;
 715     gmx_molblock_t *molb;
 716     gmx_moltype_t  *molt;
 717
 718     if (fplog)
 719     {
 720         fprintf(fplog, "\nLinking all bonded interactions to atoms\n");
 721     }
 722
 723     /* If normal and/or settle constraints act only within charge groups,
 724      * we can store them in the reverse top and simply assign them to domains.
 725      * Otherwise we need to assign them to multiple domains and set up
 726      * the parallel version constraint algoirthm(s).
 727      */
 728
 729     dd->reverse_top = make_reverse_top(mtop, ir->efep != efepNO,
 730                                        vsite ? vsite->vsite_pbc_molt : NULL,
 731                                        !dd->bInterCGcons, !dd->bInterCGsettles,
 732                                        bBCheck, &dd->nbonded_global);
 733
 734     if (dd->reverse_top->ril_mt_tot_size >= 200000 &&
 735         mtop->mols.nr > 1 &&
 736         mtop->nmolblock == 1 && mtop->molblock[0].nmol == 1)
 737     {
 738         /* mtop comes from a pre Gromacs 4 tpr file */
 739         const char *note = "NOTE: The tpr file used for this simulation is in an old format, for less memory usage and possibly more performance create a new tpr file with an up to date version of grompp";
 740         if (fplog)
 741         {
 742             fprintf(fplog, "\n%s\n\n", note);
 743         }
 744         if (DDMASTER(dd))
 745         {
 746             fprintf(stderr, "\n%s\n\n", note);
 747         }
 748     }
 749
 750     dd->reverse_top->bExclRequired = IR_EXCL_FORCES(*ir);
 751
 752     nexcl              = 0;
 753     dd->n_intercg_excl = 0;
 754     for (mb = 0; mb < mtop->nmolblock; mb++)
 755     {
 756         molb                = &mtop->molblock[mb];
 757         molt                = &mtop->moltype[molb->type];
 758         nexcl              += molb->nmol*count_excls(&molt->cgs, &molt->excls, &nexcl_icg);
 759         dd->n_intercg_excl += molb->nmol*nexcl_icg;
 760     }
 761     if (dd->reverse_top->bExclRequired)
 762     {
 763         dd->nbonded_global += nexcl;
 764         if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0 && fplog)
 765         {
 766             fprintf(fplog, "There are %d inter charge-group exclusions,\n"
 767                     "will use an extra communication step for exclusion forces for %s\n",
 768                     dd->n_intercg_excl, eel_names[ir->coulombtype]);
 769         }
 770     }
 771
 772     if (vsite && vsite->n_intercg_vsite > 0)
 773     {
 774         if (fplog)
 775         {
 776             fprintf(fplog, "There are %d inter charge-group virtual sites,\n"
 777                     "will an extra communication step for selected coordinates and forces\n",
 778                     vsite->n_intercg_vsite);
 779         }
 780         init_domdec_vsites(dd, vsite->n_intercg_vsite);
 781     }
 782
 783     if (dd->bInterCGcons || dd->bInterCGsettles)
 784     {
 785         init_domdec_constraints(dd, mtop);
 786     }
 787     if (fplog)
 788     {
 789         fprintf(fplog, "\n");
 790     }
 791 }
 792
 793 /* This routine is very similar to add_ifunc, but vsites interactions
 794  * have more work to do than other kinds of interactions, and the
 795  * complex way nral (and thus vector contents) depends on ftype
 796  * confuses static analysis tools unless we fuse the vsite
 797  * atom-indexing organization code with the ifunc-adding code, so that
 798  * they can see that nral is the same value. */
 799 static gmx_inline void
 800 add_ifunc_for_vsites(t_iatom *tiatoms, gmx_ga2la_t ga2la,
 801                      int nral, gmx_bool bHomeA,
 802                      int a, int a_gl, int a_mol,
 803                      const t_iatom *iatoms,
 804                      t_ilist *il)
 805 {
 806     t_iatom *liatoms;
 807
 808     if (il->nr+1+nral > il->nalloc)
 809     {
 810         il->nalloc = over_alloc_large(il->nr+1+nral);
 811         srenew(il->iatoms, il->nalloc);
 812     }
 813     liatoms = il->iatoms + il->nr;
 814     il->nr += 1 + nral;
 815
 816     /* Copy the type */
 817     tiatoms[0] = iatoms[0];
 818
 819     if (bHomeA)
 820     {
 821         /* We know the local index of the first atom */
 822         tiatoms[1] = a;
 823     }
 824     else
 825     {
 826         /* Convert later in make_local_vsites */
 827         tiatoms[1] = -a_gl - 1;
 828     }
 829
 830     for (int k = 2; k < 1+nral; k++)
 831     {
 832         int ak_gl = a_gl + iatoms[k] - a_mol;
 833         if (!ga2la_get_home(ga2la, ak_gl, &tiatoms[k]))
 834         {
 835             /* Copy the global index, convert later in make_local_vsites */
 836             tiatoms[k] = -(ak_gl + 1);
 837         }
 838         // Note that ga2la_get_home always sets the third parameter if
 839         // it returns TRUE
 840     }
 841     for (int k = 0; k < 1+nral; k++)
 842     {
 843         liatoms[k] = tiatoms[k];
 844     }
 845 }
 846
 847 static gmx_inline void add_ifunc(int nral, t_iatom *tiatoms, t_ilist *il)
 848 {
 849     t_iatom *liatoms;
 850     int      k;
 851
 852     if (il->nr+1+nral > il->nalloc)
 853     {
 854         il->nalloc = over_alloc_large(il->nr+1+nral);
 855         srenew(il->iatoms, il->nalloc);
 856     }
 857     liatoms = il->iatoms + il->nr;
 858     for (k = 0; k <= nral; k++)
 859     {
 860         liatoms[k] = tiatoms[k];
 861     }
 862     il->nr += 1 + nral;
 863 }
 864
 865 static void add_posres(int mol, int a_mol, const gmx_molblock_t *molb,
 866                        t_iatom *iatoms, const t_iparams *ip_in,
 867                        t_idef *idef)
 868 {
 869     int        n, a_molb;
 870     t_iparams *ip;
 871
 872     /* This position restraint has not been added yet,
 873      * so it's index is the current number of position restraints.
 874      */
 875     n = idef->il[F_POSRES].nr/2;
 876     if (n+1 > idef->iparams_posres_nalloc)
 877     {
 878         idef->iparams_posres_nalloc = over_alloc_dd(n+1);
 879         srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
 880     }
 881     ip = &idef->iparams_posres[n];
 882     /* Copy the force constants */
 883     *ip = ip_in[iatoms[0]];
 884
 885     /* Get the position restraint coordinates from the molblock */
 886     a_molb = mol*molb->natoms_mol + a_mol;
 887     if (a_molb >= molb->nposres_xA)
 888     {
 889         gmx_incons("Not enough position restraint coordinates");
 890     }
 891     ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 892     ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
 893     ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
 894     if (molb->nposres_xB > 0)
 895     {
 896         ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
 897         ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
 898         ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
 899     }
 900     else
 901     {
 902         ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
 903         ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
 904         ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
 905     }
 906     /* Set the parameter index for idef->iparams_posre */
 907     iatoms[0] = n;
 908 }
 909
 910 static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb,
 911                          t_iatom *iatoms, const t_iparams *ip_in,
 912                          t_idef *idef)
 913 {
 914     int        n, a_molb;
 915     t_iparams *ip;
 916
 917     /* This flat-bottom position restraint has not been added yet,
 918      * so it's index is the current number of position restraints.
 919      */
 920     n = idef->il[F_FBPOSRES].nr/2;
 921     if (n+1 > idef->iparams_fbposres_nalloc)
 922     {
 923         idef->iparams_fbposres_nalloc = over_alloc_dd(n+1);
 924         srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
 925     }
 926     ip = &idef->iparams_fbposres[n];
 927     /* Copy the force constants */
 928     *ip = ip_in[iatoms[0]];
 929
 930     /* Get the position restriant coordinats from the molblock */
 931     a_molb = mol*molb->natoms_mol + a_mol;
 932     if (a_molb >= molb->nposres_xA)
 933     {
 934         gmx_incons("Not enough position restraint coordinates");
 935     }
 936     /* Take reference positions from A position of normal posres */
 937     ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 938     ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 939     ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 940
 941     /* Note: no B-type for flat-bottom posres */
 942
 943     /* Set the parameter index for idef->iparams_posre */
 944     iatoms[0] = n;
 945 }
 946
 947 static void add_vsite(gmx_ga2la_t ga2la, int *index, int *rtil,
 948                       int ftype, int nral,
 949                       gmx_bool bHomeA, int a, int a_gl, int a_mol,
 950                       t_iatom *iatoms,
 951                       t_idef *idef, int **vsite_pbc, int *vsite_pbc_nalloc)
 952 {
 953     int     k, vsi, pbc_a_mol;
 954     t_iatom tiatoms[1+MAXATOMLIST];
 955     int     j, ftype_r, nral_r;
 956
 957     /* Add this interaction to the local topology */
 958     add_ifunc_for_vsites(tiatoms, ga2la, nral, bHomeA, a, a_gl, a_mol, iatoms, &idef->il[ftype]);
 959
 960     if (vsite_pbc)
 961     {
 962         vsi = idef->il[ftype].nr/(1+nral) - 1;
 963         if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
 964         {
 965             vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
 966             srenew(vsite_pbc[ftype-F_VSITE2], vsite_pbc_nalloc[ftype-F_VSITE2]);
 967         }
 968         if (bHomeA)
 969         {
 970             pbc_a_mol = iatoms[1+nral+1];
 971             if (pbc_a_mol < 0)
 972             {
 973                 /* The pbc flag is one of the following two options:
 974                  * -2: vsite and all constructing atoms are within the same cg, no pbc
 975                  * -1: vsite and its first constructing atom are in the same cg, do pbc
 976                  */
 977                 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
 978             }
 979             else
 980             {
 981                 /* Set the pbc atom for this vsite so we can make its pbc
 982                  * identical to the rest of the atoms in its charge group.
 983                  * Since the order of the atoms does not change within a charge
 984                  * group, we do not need the global to local atom index.
 985                  */
 986                 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
 987             }
 988         }
 989         else
 990         {
 991             /* This vsite is non-home (required for recursion),
 992              * and therefore there is no charge group to match pbc with.
 993              * But we always turn on full_pbc to assure that higher order
 994              * recursion works correctly.
 995              */
 996             vsite_pbc[ftype-F_VSITE2][vsi] = -1;
 997         }
 998     }
 999
1000     if (iatoms[1+nral])
1001     {
1002         /* Check for recursion */
1003         for (k = 2; k < 1+nral; k++)
1004         {
1005             if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
1006             {
1007                 /* This construction atoms is a vsite and not a home atom */
1008                 if (gmx_debug_at)
1009                 {
1010                     fprintf(debug, "Constructing atom %d of vsite atom %d is a vsite and non-home\n", iatoms[k]+1, a_mol+1);
1011                 }
1012                 /* Find the vsite construction */
1013
1014                 /* Check all interactions assigned to this atom */
1015                 j = index[iatoms[k]];
1016                 while (j < index[iatoms[k]+1])
1017                 {
1018                     ftype_r = rtil[j++];
1019                     nral_r  = NRAL(ftype_r);
1020                     if (interaction_function[ftype_r].flags & IF_VSITE)
1021                     {
1022                         /* Add this vsite (recursion) */
1023                         add_vsite(ga2la, index, rtil, ftype_r, nral_r,
1024                                   FALSE, -1, a_gl+iatoms[k]-iatoms[1], iatoms[k],
1025                                   rtil+j, idef, vsite_pbc, vsite_pbc_nalloc);
1026                         j += 1 + nral_r + 2;
1027                     }
1028                     else
1029                     {
1030                         j += 1 + nral_r;
1031                     }
1032                 }
1033             }
1034         }
1035     }
1036 }
1037
1038 static void make_la2lc(gmx_domdec_t *dd)
1039 {
1040     int *cgindex, *la2lc, cg, a;
1041
1042     cgindex = dd->cgindex;
1043
1044     if (dd->nat_tot > dd->la2lc_nalloc)
1045     {
1046         dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
1047         snew(dd->la2lc, dd->la2lc_nalloc);
1048     }
1049     la2lc = dd->la2lc;
1050
1051     /* Make the local atom to local cg index */
1052     for (cg = 0; cg < dd->ncg_tot; cg++)
1053     {
1054         for (a = cgindex[cg]; a < cgindex[cg+1]; a++)
1055         {
1056             la2lc[a] = cg;
1057         }
1058     }
1059 }
1060
1061 static real dd_dist2(t_pbc *pbc_null, rvec *cg_cm, const int *la2lc, int i, int j)
1062 {
1063     rvec dx;
1064
1065     if (pbc_null)
1066     {
1067         pbc_dx_aiuc(pbc_null, cg_cm[la2lc[i]], cg_cm[la2lc[j]], dx);
1068     }
1069     else
1070     {
1071         rvec_sub(cg_cm[la2lc[i]], cg_cm[la2lc[j]], dx);
1072     }
1073
1074     return norm2(dx);
1075 }
1076
1077 /* Append the nsrc t_blocka block structures in src to *dest */
1078 static void combine_blocka(t_blocka *dest, const t_blocka *src, int nsrc)
1079 {
1080     int ni, na, s, i;
1081
1082     ni = src[nsrc-1].nr;
1083     na = 0;
1084     for (s = 0; s < nsrc; s++)
1085     {
1086         na += src[s].nra;
1087     }
1088     if (ni + 1 > dest->nalloc_index)
1089     {
1090         dest->nalloc_index = over_alloc_large(ni+1);
1091         srenew(dest->index, dest->nalloc_index);
1092     }
1093     if (dest->nra + na > dest->nalloc_a)
1094     {
1095         dest->nalloc_a = over_alloc_large(dest->nra+na);
1096         srenew(dest->a, dest->nalloc_a);
1097     }
1098     for (s = 0; s < nsrc; s++)
1099     {
1100         for (i = dest->nr+1; i < src[s].nr+1; i++)
1101         {
1102             dest->index[i] = dest->nra + src[s].index[i];
1103         }
1104         for (i = 0; i < src[s].nra; i++)
1105         {
1106             dest->a[dest->nra+i] = src[s].a[i];
1107         }
1108         dest->nr   = src[s].nr;
1109         dest->nra += src[s].nra;
1110     }
1111 }
1112
1113 /* Append the nsrc t_idef structures in src to *dest,
1114  * virtual sites need special attention, as pbc info differs per vsite.
1115  */
1116 static void combine_idef(t_idef *dest, const t_idef *src, int nsrc,
1117                          gmx_vsite_t *vsite, int ***vsite_pbc_t)
1118 {
1119     int            ftype, n, s, i;
1120     t_ilist       *ild;
1121     const t_ilist *ils;
1122     gmx_bool       vpbc;
1123     int            nral1 = 0, ftv = 0;
1124
1125     for (ftype = 0; ftype < F_NRE; ftype++)
1126     {
1127         n = 0;
1128         for (s = 0; s < nsrc; s++)
1129         {
1130             n += src[s].il[ftype].nr;
1131         }
1132         if (n > 0)
1133         {
1134             ild = &dest->il[ftype];
1135
1136             if (ild->nr + n > ild->nalloc)
1137             {
1138                 ild->nalloc = over_alloc_large(ild->nr+n);
1139                 srenew(ild->iatoms, ild->nalloc);
1140             }
1141
1142             vpbc = ((interaction_function[ftype].flags & IF_VSITE) &&
1143                     vsite->vsite_pbc_loc != NULL);
1144             if (vpbc)
1145             {
1146                 nral1 = 1 + NRAL(ftype);
1147                 ftv   = ftype - F_VSITE2;
1148                 if ((ild->nr + n)/nral1 > vsite->vsite_pbc_loc_nalloc[ftv])
1149                 {
1150                     vsite->vsite_pbc_loc_nalloc[ftv] =
1151                         over_alloc_large((ild->nr + n)/nral1);
1152                     srenew(vsite->vsite_pbc_loc[ftv],
1153                            vsite->vsite_pbc_loc_nalloc[ftv]);
1154                 }
1155             }
1156
1157             for (s = 0; s < nsrc; s++)
1158             {
1159                 ils = &src[s].il[ftype];
1160                 for (i = 0; i < ils->nr; i++)
1161                 {
1162                     ild->iatoms[ild->nr+i] = ils->iatoms[i];
1163                 }
1164                 if (vpbc)
1165                 {
1166                     for (i = 0; i < ils->nr; i += nral1)
1167                     {
1168                         vsite->vsite_pbc_loc[ftv][(ild->nr+i)/nral1] =
1169                             vsite_pbc_t[s][ftv][i/nral1];
1170                     }
1171                 }
1172
1173                 ild->nr += ils->nr;
1174             }
1175         }
1176     }
1177
1178     /* Position restraints need an additional treatment */
1179     if (dest->il[F_POSRES].nr > 0)
1180     {
1181         n = dest->il[F_POSRES].nr/2;
1182         if (n > dest->iparams_posres_nalloc)
1183         {
1184             dest->iparams_posres_nalloc = over_alloc_large(n);
1185             srenew(dest->iparams_posres, dest->iparams_posres_nalloc);
1186         }
1187         /* Set n to the number of original position restraints in dest */
1188         for (s = 0; s < nsrc; s++)
1189         {
1190             n -= src[s].il[F_POSRES].nr/2;
1191         }
1192         for (s = 0; s < nsrc; s++)
1193         {
1194             for (i = 0; i < src[s].il[F_POSRES].nr/2; i++)
1195             {
1196                 /* Correct the index into iparams_posres */
1197                 dest->il[F_POSRES].iatoms[n*2] = n;
1198                 /* Copy the position restraint force parameters */
1199                 dest->iparams_posres[n] = src[s].iparams_posres[i];
1200                 n++;
1201             }
1202         }
1203     }
1204 }
1205
1206 /* This function looks up and assigns bonded interactions for zone iz.
1207  * With thread parallelizing each thread acts on a different atom range:
1208  * at_start to at_end.
1209  */
1210 static int make_bondeds_zone(gmx_domdec_t *dd,
1211                              const gmx_domdec_zones_t *zones,
1212                              const gmx_molblock_t *molb,
1213                              gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
1214                              real rc2,
1215                              int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
1216                              const t_iparams *ip_in,
1217                              t_idef *idef,
1218                              int **vsite_pbc,
1219                              int *vsite_pbc_nalloc,
1220                              int iz, int nzone,
1221                              int at_start, int at_end)
1222 {
1223     int                           i, i_gl, mb, mt, mol, i_mol, j, ftype, nral, d, k;
1224     int                          *index, *rtil;
1225     t_iatom                      *iatoms, tiatoms[1+MAXATOMLIST];
1226     gmx_bool                      bBCheck, bUse, bLocal;
1227     ivec                          k_zero, k_plus;
1228     gmx_ga2la_t                   ga2la;
1229     int                           a_loc;
1230     int                           kz;
1231     int                           nizone;
1232     const gmx_domdec_ns_ranges_t *izone;
1233     gmx_reverse_top_t            *rt;
1234     int                           nbonded_local;
1235
1236     nizone = zones->nizone;
1237     izone  = zones->izone;
1238
1239     rt = dd->reverse_top;
1240
1241     bBCheck = rt->bBCheck;
1242
1243     nbonded_local = 0;
1244
1245     ga2la = dd->ga2la;
1246
1247     for (i = at_start; i < at_end; i++)
1248     {
1249         /* Get the global atom number */
1250         i_gl = dd->gatindex[i];
1251         global_atomnr_to_moltype_ind(rt, i_gl, &mb, &mt, &mol, &i_mol);
1252         /* Check all interactions assigned to this atom */
1253         index = rt->ril_mt[mt].index;
1254         rtil  = rt->ril_mt[mt].il;
1255         j     = index[i_mol];
1256         while (j < index[i_mol+1])
1257         {
1258             ftype  = rtil[j++];
1259             iatoms = rtil + j;
1260             nral   = NRAL(ftype);
1261             if (ftype == F_SETTLE)
1262             {
1263                 /* Settles are only in the reverse top when they
1264                  * operate within a charge group. So we can assign
1265                  * them without checks. We do this only for performance
1266                  * reasons; it could be handled by the code below.
1267                  */
1268                 if (iz == 0)
1269                 {
1270                     /* Home zone: add this settle to the local topology */
1271                     tiatoms[0] = iatoms[0];
1272                     tiatoms[1] = i;
1273                     tiatoms[2] = i + iatoms[2] - iatoms[1];
1274                     tiatoms[3] = i + iatoms[3] - iatoms[1];
1275                     add_ifunc(nral, tiatoms, &idef->il[ftype]);
1276                     nbonded_local++;
1277                 }
1278                 j += 1 + nral;
1279             }
1280             else if (interaction_function[ftype].flags & IF_VSITE)
1281             {
1282                 /* The vsite construction goes where the vsite itself is */
1283                 if (iz == 0)
1284                 {
1285                     add_vsite(dd->ga2la, index, rtil, ftype, nral,
1286                               TRUE, i, i_gl, i_mol,
1287                               iatoms, idef, vsite_pbc, vsite_pbc_nalloc);
1288                 }
1289                 j += 1 + nral + 2;
1290             }
1291             else
1292             {
1293                 /* Copy the type */
1294                 tiatoms[0] = iatoms[0];
1295
1296                 if (nral == 1)
1297                 {
1298                     /* Assign single-body interactions to the home zone */
1299                     if (iz == 0)
1300                     {
1301                         bUse       = TRUE;
1302                         tiatoms[1] = i;
1303                         if (ftype == F_POSRES)
1304                         {
1305                             add_posres(mol, i_mol, &molb[mb], tiatoms, ip_in,
1306                                        idef);
1307                         }
1308                         else if (ftype == F_FBPOSRES)
1309                         {
1310                             add_fbposres(mol, i_mol, &molb[mb], tiatoms, ip_in,
1311                                          idef);
1312                         }
1313                     }
1314                     else
1315                     {
1316                         bUse = FALSE;
1317                     }
1318                 }
1319                 else if (nral == 2)
1320                 {
1321                     /* This is a two-body interaction, we can assign
1322                      * analogous to the non-bonded assignments.
1323                      */
1324                     if (!ga2la_get(ga2la, i_gl+iatoms[2]-i_mol, &a_loc, &kz))
1325                     {
1326                         bUse = FALSE;
1327                     }
1328                     else
1329                     {
1330                         if (kz >= nzone)
1331                         {
1332                             kz -= nzone;
1333                         }
1334                         /* Check zone interaction assignments */
1335                         bUse = ((iz < nizone && iz <= kz &&
1336                                  izone[iz].j0 <= kz && kz < izone[iz].j1) ||
1337                                 (kz < nizone &&iz >  kz &&
1338                                  izone[kz].j0 <= iz && iz < izone[kz].j1));
1339                         if (bUse)
1340                         {
1341                             tiatoms[1] = i;
1342                             tiatoms[2] = a_loc;
1343                             /* If necessary check the cgcm distance */
1344                             if (bRCheck2B &&
1345                                 dd_dist2(pbc_null, cg_cm, la2lc,
1346                                          tiatoms[1], tiatoms[2]) >= rc2)
1347                             {
1348                                 bUse = FALSE;
1349                             }
1350                         }
1351                     }
1352                 }
1353                 else
1354                 {
1355                     /* Assign this multi-body bonded interaction to
1356                      * the local node if we have all the atoms involved
1357                      * (local or communicated) and the minimum zone shift
1358                      * in each dimension is zero, for dimensions
1359                      * with 2 DD cells an extra check may be necessary.
1360                      */
1361                     bUse = TRUE;
1362                     clear_ivec(k_zero);
1363                     clear_ivec(k_plus);
1364                     for (k = 1; k <= nral && bUse; k++)
1365                     {
1366                         bLocal = ga2la_get(ga2la, i_gl+iatoms[k]-i_mol,
1367                                            &a_loc, &kz);
1368                         if (!bLocal || kz >= zones->n)
1369                         {
1370                             /* We do not have this atom of this interaction
1371                              * locally, or it comes from more than one cell
1372                              * away.
1373                              */
1374                             bUse = FALSE;
1375                         }
1376                         else
1377                         {
1378                             tiatoms[k] = a_loc;
1379                             for (d = 0; d < DIM; d++)
1380                             {
1381                                 if (zones->shift[kz][d] == 0)
1382                                 {
1383                                     k_zero[d] = k;
1384                                 }
1385                                 else
1386                                 {
1387                                     k_plus[d] = k;
1388                                 }
1389                             }
1390                         }
1391                     }
1392                     bUse = (bUse &&
1393                             k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1394                     if (bRCheckMB)
1395                     {
1396                         for (d = 0; (d < DIM && bUse); d++)
1397                         {
1398                             /* Check if the cg_cm distance falls within
1399                              * the cut-off to avoid possible multiple
1400                              * assignments of bonded interactions.
1401                              */
1402                             if (rcheck[d] &&
1403                                 k_plus[d] &&
1404                                 dd_dist2(pbc_null, cg_cm, la2lc,
1405                                          tiatoms[k_zero[d]], tiatoms[k_plus[d]]) >= rc2)
1406                             {
1407                                 bUse = FALSE;
1408                             }
1409                         }
1410                     }
1411                 }
1412                 if (bUse)
1413                 {
1414                     /* Add this interaction to the local topology */
1415                     add_ifunc(nral, tiatoms, &idef->il[ftype]);
1416                     /* Sum so we can check in global_stat
1417                      * if we have everything.
1418                      */
1419                     if (bBCheck ||
1420                         !(interaction_function[ftype].flags & IF_LIMZERO))
1421                     {
1422                         nbonded_local++;
1423                     }
1424                 }
1425                 j += 1 + nral;
1426             }
1427         }
1428     }
1429
1430     return nbonded_local;
1431 }
1432
1433 static void set_no_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1434                                    int iz, t_blocka *lexcls)
1435 {
1436     int  a0, a1, a;
1437
1438     a0 = dd->cgindex[zones->cg_range[iz]];
1439     a1 = dd->cgindex[zones->cg_range[iz+1]];
1440
1441     for (a = a0+1; a < a1+1; a++)
1442     {
1443         lexcls->index[a] = lexcls->nra;
1444     }
1445 }
1446
1447 static int make_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1448                                 const gmx_moltype_t *moltype,
1449                                 gmx_bool bRCheck, real rc2,
1450                                 int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
1451                                 const int *cginfo,
1452                                 t_blocka *lexcls,
1453                                 int iz,
1454                                 int cg_start, int cg_end)
1455 {
1456     int             n, count, jla0, jla1, jla;
1457     int             cg, la0, la1, la, a_gl, mb, mt, mol, a_mol, j, aj_mol;
1458     const t_blocka *excls;
1459     gmx_ga2la_t     ga2la;
1460     int             cell;
1461
1462     ga2la = dd->ga2la;
1463
1464     jla0 = dd->cgindex[zones->izone[iz].jcg0];
1465     jla1 = dd->cgindex[zones->izone[iz].jcg1];
1466
1467     /* We set the end index, but note that we might not start at zero here */
1468     lexcls->nr = dd->cgindex[cg_end];
1469
1470     n     = lexcls->nra;
1471     count = 0;
1472     for (cg = cg_start; cg < cg_end; cg++)
1473     {
1474         /* Here we assume the number of exclusions in one charge group
1475          * is never larger than 1000.
1476          */
1477         if (n+1000 > lexcls->nalloc_a)
1478         {
1479             lexcls->nalloc_a = over_alloc_large(n+1000);
1480             srenew(lexcls->a, lexcls->nalloc_a);
1481         }
1482         la0 = dd->cgindex[cg];
1483         la1 = dd->cgindex[cg+1];
1484         if (GET_CGINFO_EXCL_INTER(cginfo[cg]) ||
1485             !GET_CGINFO_EXCL_INTRA(cginfo[cg]))
1486         {
1487             /* Copy the exclusions from the global top */
1488             for (la = la0; la < la1; la++)
1489             {
1490                 lexcls->index[la] = n;
1491                 a_gl              = dd->gatindex[la];
1492                 global_atomnr_to_moltype_ind(dd->reverse_top, a_gl, &mb, &mt, &mol, &a_mol);
1493                 excls = &moltype[mt].excls;
1494                 for (j = excls->index[a_mol]; j < excls->index[a_mol+1]; j++)
1495                 {
1496                     aj_mol = excls->a[j];
1497                     /* This computation of jla is only correct intra-cg */
1498                     jla = la + aj_mol - a_mol;
1499                     if (jla >= la0 && jla < la1)
1500                     {
1501                         /* This is an intra-cg exclusion. We can skip
1502                          *  the global indexing and distance checking.
1503                          */
1504                         /* Intra-cg exclusions are only required
1505                          * for the home zone.
1506                          */
1507                         if (iz == 0)
1508                         {
1509                             lexcls->a[n++] = jla;
1510                             /* Check to avoid double counts */
1511                             if (jla > la)
1512                             {
1513                                 count++;
1514                             }
1515                         }
1516                     }
1517                     else
1518                     {
1519                         /* This is a inter-cg exclusion */
1520                         /* Since exclusions are pair interactions,
1521                          * just like non-bonded interactions,
1522                          * they can be assigned properly up
1523                          * to the DD cutoff (not cutoff_min as
1524                          * for the other bonded interactions).
1525                          */
1526                         if (ga2la_get(ga2la, a_gl+aj_mol-a_mol, &jla, &cell))
1527                         {
1528                             if (iz == 0 && cell == 0)
1529                             {
1530                                 lexcls->a[n++] = jla;
1531                                 /* Check to avoid double counts */
1532                                 if (jla > la)
1533                                 {
1534                                     count++;
1535                                 }
1536                             }
1537                             else if (jla >= jla0 && jla < jla1 &&
1538                                      (!bRCheck ||
1539                                       dd_dist2(pbc_null, cg_cm, la2lc, la, jla) < rc2))
1540                             {
1541                                 /* jla > la, since jla0 > la */
1542                                 lexcls->a[n++] = jla;
1543                                 count++;
1544                             }
1545                         }
1546                     }
1547                 }
1548             }
1549         }
1550         else
1551         {
1552             /* There are no inter-cg excls and this cg is self-excluded.
1553              * These exclusions are only required for zone 0,
1554              * since other zones do not see themselves.
1555              */
1556             if (iz == 0)
1557             {
1558                 for (la = la0; la < la1; la++)
1559                 {
1560                     lexcls->index[la] = n;
1561                     for (j = la0; j < la1; j++)
1562                     {
1563                         lexcls->a[n++] = j;
1564                     }
1565                 }
1566                 count += ((la1 - la0)*(la1 - la0 - 1))/2;
1567             }
1568             else
1569             {
1570                 /* We don't need exclusions for this cg */
1571                 for (la = la0; la < la1; la++)
1572                 {
1573                     lexcls->index[la] = n;
1574                 }
1575             }
1576         }
1577     }
1578
1579     lexcls->index[lexcls->nr] = n;
1580     lexcls->nra               = n;
1581
1582     return count;
1583 }
1584
1585 static void check_alloc_index(t_blocka *ba, int nindex_max)
1586 {
1587     if (nindex_max+1 > ba->nalloc_index)
1588     {
1589         ba->nalloc_index = over_alloc_dd(nindex_max+1);
1590         srenew(ba->index, ba->nalloc_index);
1591     }
1592 }
1593
1594 static void check_exclusions_alloc(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1595                                    t_blocka *lexcls)
1596 {
1597     int nr;
1598     int thread;
1599
1600     nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1601
1602     check_alloc_index(lexcls, nr);
1603
1604     for (thread = 1; thread < dd->reverse_top->nthread; thread++)
1605     {
1606         check_alloc_index(&dd->reverse_top->excl_thread[thread], nr);
1607     }
1608 }
1609
1610 static void finish_local_exclusions(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1611                                     t_blocka *lexcls)
1612 {
1613     int la0, la;
1614
1615     lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1616
1617     if (dd->n_intercg_excl == 0)
1618     {
1619         /* There are no exclusions involving non-home charge groups,
1620          * but we need to set the indices for neighborsearching.
1621          */
1622         la0 = dd->cgindex[zones->izone[0].cg1];
1623         for (la = la0; la < lexcls->nr; la++)
1624         {
1625             lexcls->index[la] = lexcls->nra;
1626         }
1627
1628         /* nr is only used to loop over the exclusions for Ewald and RF,
1629          * so we can set it to the number of home atoms for efficiency.
1630          */
1631         lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1632     }
1633 }
1634
1635 static void clear_idef(t_idef *idef)
1636 {
1637     int  ftype;
1638
1639     /* Clear the counts */
1640     for (ftype = 0; ftype < F_NRE; ftype++)
1641     {
1642         idef->il[ftype].nr = 0;
1643     }
1644 }
1645
1646 static int make_local_bondeds_excls(gmx_domdec_t *dd,
1647                                     gmx_domdec_zones_t *zones,
1648                                     const gmx_mtop_t *mtop,
1649                                     const int *cginfo,
1650                                     gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
1651                                     real rc,
1652                                     int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
1653                                     t_idef *idef, gmx_vsite_t *vsite,
1654                                     t_blocka *lexcls, int *excl_count)
1655 {
1656     int                nzone_bondeds, nzone_excl;
1657     int                iz, cg0, cg1;
1658     real               rc2;
1659     int                nbonded_local;
1660     int                thread;
1661     gmx_reverse_top_t *rt;
1662
1663     if (dd->reverse_top->bMultiCGmols)
1664     {
1665         nzone_bondeds = zones->n;
1666     }
1667     else
1668     {
1669         /* Only single charge group molecules, so interactions don't
1670          * cross zone boundaries and we only need to assign in the home zone.
1671          */
1672         nzone_bondeds = 1;
1673     }
1674
1675     if (dd->n_intercg_excl > 0)
1676     {
1677         /* We only use exclusions from i-zones to i- and j-zones */
1678         nzone_excl = zones->nizone;
1679     }
1680     else
1681     {
1682         /* There are no inter-cg exclusions and only zone 0 sees itself */
1683         nzone_excl = 1;
1684     }
1685
1686     check_exclusions_alloc(dd, zones, lexcls);
1687
1688     rt = dd->reverse_top;
1689
1690     rc2 = rc*rc;
1691
1692     /* Clear the counts */
1693     clear_idef(idef);
1694     nbonded_local = 0;
1695
1696     lexcls->nr    = 0;
1697     lexcls->nra   = 0;
1698     *excl_count   = 0;
1699
1700     for (iz = 0; iz < nzone_bondeds; iz++)
1701     {
1702         cg0 = zones->cg_range[iz];
1703         cg1 = zones->cg_range[iz+1];
1704
1705 #pragma omp parallel for num_threads(rt->nthread) schedule(static)
1706         for (thread = 0; thread < rt->nthread; thread++)
1707         {
1708             int       cg0t, cg1t;
1709             t_idef   *idef_t;
1710             int     **vsite_pbc;
1711             int      *vsite_pbc_nalloc;
1712             t_blocka *excl_t;
1713
1714             cg0t = cg0 + ((cg1 - cg0)* thread   )/rt->nthread;
1715             cg1t = cg0 + ((cg1 - cg0)*(thread+1))/rt->nthread;
1716
1717             if (thread == 0)
1718             {
1719                 idef_t = idef;
1720             }
1721             else
1722             {
1723                 idef_t = &rt->idef_thread[thread];
1724                 clear_idef(idef_t);
1725             }
1726
1727             if (vsite && vsite->bHaveChargeGroups && vsite->n_intercg_vsite > 0)
1728             {
1729                 if (thread == 0)
1730                 {
1731                     vsite_pbc        = vsite->vsite_pbc_loc;
1732                     vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1733                 }
1734                 else
1735                 {
1736                     vsite_pbc        = rt->vsite_pbc[thread];
1737                     vsite_pbc_nalloc = rt->vsite_pbc_nalloc[thread];
1738                 }
1739             }
1740             else
1741             {
1742                 vsite_pbc        = NULL;
1743                 vsite_pbc_nalloc = NULL;
1744             }
1745
1746             rt->nbonded_thread[thread] =
1747                 make_bondeds_zone(dd, zones,
1748                                   mtop->molblock,
1749                                   bRCheckMB, rcheck, bRCheck2B, rc2,
1750                                   la2lc, pbc_null, cg_cm, idef->iparams,
1751                                   idef_t,
1752                                   vsite_pbc, vsite_pbc_nalloc,
1753                                   iz, zones->n,
1754                                   dd->cgindex[cg0t], dd->cgindex[cg1t]);
1755
1756             if (iz < nzone_excl)
1757             {
1758                 if (thread == 0)
1759                 {
1760                     excl_t = lexcls;
1761                 }
1762                 else
1763                 {
1764                     excl_t      = &rt->excl_thread[thread];
1765                     excl_t->nr  = 0;
1766                     excl_t->nra = 0;
1767                 }
1768
1769                 rt->excl_count_thread[thread] =
1770                     make_exclusions_zone(dd, zones,
1771                                          mtop->moltype, bRCheck2B, rc2,
1772                                          la2lc, pbc_null, cg_cm, cginfo,
1773                                          excl_t,
1774                                          iz,
1775                                          cg0t, cg1t);
1776             }
1777         }
1778
1779         if (rt->nthread > 1)
1780         {
1781             combine_idef(idef, rt->idef_thread+1, rt->nthread-1,
1782                          vsite, rt->vsite_pbc+1);
1783         }
1784
1785         for (thread = 0; thread < rt->nthread; thread++)
1786         {
1787             nbonded_local += rt->nbonded_thread[thread];
1788         }
1789
1790         if (iz < nzone_excl)
1791         {
1792             if (rt->nthread > 1)
1793             {
1794                 combine_blocka(lexcls, rt->excl_thread+1, rt->nthread-1);
1795             }
1796
1797             for (thread = 0; thread < rt->nthread; thread++)
1798             {
1799                 *excl_count += rt->excl_count_thread[thread];
1800             }
1801         }
1802     }
1803
1804     /* Some zones might not have exclusions, but some code still needs to
1805      * loop over the index, so we set the indices here.
1806      */
1807     for (iz = nzone_excl; iz < zones->nizone; iz++)
1808     {
1809         set_no_exclusions_zone(dd, zones, iz, lexcls);
1810     }
1811
1812     finish_local_exclusions(dd, zones, lexcls);
1813     if (debug)
1814     {
1815         fprintf(debug, "We have %d exclusions, check count %d\n",
1816                 lexcls->nra, *excl_count);
1817     }
1818
1819     return nbonded_local;
1820 }
1821
1822 void dd_make_local_cgs(gmx_domdec_t *dd, t_block *lcgs)
1823 {
1824     lcgs->nr    = dd->ncg_tot;
1825     lcgs->index = dd->cgindex;
1826 }
1827
1828 void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
1829                        int npbcdim, matrix box,
1830                        rvec cellsize_min, ivec npulse,
1831                        t_forcerec *fr,
1832                        rvec *cgcm_or_x,
1833                        gmx_vsite_t *vsite,
1834                        gmx_mtop_t *mtop, gmx_localtop_t *ltop)
1835 {
1836     gmx_bool bRCheckMB, bRCheck2B;
1837     real     rc = -1;
1838     ivec     rcheck;
1839     int      d, nexcl;
1840     t_pbc    pbc, *pbc_null = NULL;
1841
1842     if (debug)
1843     {
1844         fprintf(debug, "Making local topology\n");
1845     }
1846
1847     dd_make_local_cgs(dd, &ltop->cgs);
1848
1849     bRCheckMB   = FALSE;
1850     bRCheck2B   = FALSE;
1851
1852     if (dd->reverse_top->bMultiCGmols)
1853     {
1854         /* We need to check to which cell bondeds should be assigned */
1855         rc = dd_cutoff_twobody(dd);
1856         if (debug)
1857         {
1858             fprintf(debug, "Two-body bonded cut-off distance is %g\n", rc);
1859         }
1860
1861         /* Should we check cg_cm distances when assigning bonded interactions? */
1862         for (d = 0; d < DIM; d++)
1863         {
1864             rcheck[d] = FALSE;
1865             /* Only need to check for dimensions where the part of the box
1866              * that is not communicated is smaller than the cut-off.
1867              */
1868             if (d < npbcdim && dd->nc[d] > 1 &&
1869                 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
1870             {
1871                 if (dd->nc[d] == 2)
1872                 {
1873                     rcheck[d] = TRUE;
1874                     bRCheckMB = TRUE;
1875                 }
1876                 /* Check for interactions between two atoms,
1877                  * where we can allow interactions up to the cut-off,
1878                  * instead of up to the smallest cell dimension.
1879                  */
1880                 bRCheck2B = TRUE;
1881             }
1882             if (debug)
1883             {
1884                 fprintf(debug,
1885                         "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1886                         d, cellsize_min[d], d, rcheck[d], bRCheck2B);
1887             }
1888         }
1889         if (bRCheckMB || bRCheck2B)
1890         {
1891             make_la2lc(dd);
1892             if (fr->bMolPBC)
1893             {
1894                 set_pbc_dd(&pbc, fr->ePBC, dd, TRUE, box);
1895                 pbc_null = &pbc;
1896             }
1897             else
1898             {
1899                 pbc_null = NULL;
1900             }
1901         }
1902     }
1903
1904     dd->nbonded_local =
1905         make_local_bondeds_excls(dd, zones, mtop, fr->cginfo,
1906                                  bRCheckMB, rcheck, bRCheck2B, rc,
1907                                  dd->la2lc,
1908                                  pbc_null, cgcm_or_x,
1909                                  &ltop->idef, vsite,
1910                                  &ltop->excls, &nexcl);
1911
1912     /* The ilist is not sorted yet,
1913      * we can only do this when we have the charge arrays.
1914      */
1915     ltop->idef.ilsort = ilsortUNKNOWN;
1916
1917     if (dd->reverse_top->bExclRequired)
1918     {
1919         dd->nbonded_local += nexcl;
1920
1921         forcerec_set_excl_load(fr, ltop);
1922     }
1923
1924     ltop->atomtypes  = mtop->atomtypes;
1925
1926     /* For an error message only */
1927     dd->reverse_top->err_top_global = mtop;
1928     dd->reverse_top->err_top_local  = ltop;
1929 }
1930
1931 void dd_sort_local_top(gmx_domdec_t *dd, t_mdatoms *mdatoms,
1932                        gmx_localtop_t *ltop)
1933 {
1934     if (dd->reverse_top->ilsort == ilsortNO_FE)
1935     {
1936         ltop->idef.ilsort = ilsortNO_FE;
1937     }
1938     else
1939     {
1940         gmx_sort_ilist_fe(&ltop->idef, mdatoms->chargeA, mdatoms->chargeB);
1941     }
1942 }
1943
1944 gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
1945 {
1946     gmx_localtop_t *top;
1947     int             i;
1948
1949     snew(top, 1);
1950
1951     top->idef.ntypes    = top_global->ffparams.ntypes;
1952     top->idef.atnr      = top_global->ffparams.atnr;
1953     top->idef.functype  = top_global->ffparams.functype;
1954     top->idef.iparams   = top_global->ffparams.iparams;
1955     top->idef.fudgeQQ   = top_global->ffparams.fudgeQQ;
1956     top->idef.cmap_grid = top_global->ffparams.cmap_grid;
1957
1958     for (i = 0; i < F_NRE; i++)
1959     {
1960         top->idef.il[i].iatoms = NULL;
1961         top->idef.il[i].nalloc = 0;
1962     }
1963     top->idef.ilsort   = ilsortUNKNOWN;
1964
1965     return top;
1966 }
1967
1968 void dd_init_local_state(gmx_domdec_t *dd,
1969                          t_state *state_global, t_state *state_local)
1970 {
1971     int buf[NITEM_DD_INIT_LOCAL_STATE];
1972
1973     if (DDMASTER(dd))
1974     {
1975         buf[0] = state_global->flags;
1976         buf[1] = state_global->ngtc;
1977         buf[2] = state_global->nnhpres;
1978         buf[3] = state_global->nhchainlength;
1979         buf[4] = state_global->dfhist.nlambda;
1980     }
1981     dd_bcast(dd, NITEM_DD_INIT_LOCAL_STATE*sizeof(int), buf);
1982
1983     init_state(state_local, 0, buf[1], buf[2], buf[3], buf[4]);
1984     state_local->flags = buf[0];
1985 }
1986
1987 static void check_link(t_blocka *link, int cg_gl, int cg_gl_j)
1988 {
1989     int      k;
1990     gmx_bool bFound;
1991
1992     bFound = FALSE;
1993     for (k = link->index[cg_gl]; k < link->index[cg_gl+1]; k++)
1994     {
1995         GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
1996         if (link->a[k] == cg_gl_j)
1997         {
1998             bFound = TRUE;
1999         }
2000     }
2001     if (!bFound)
2002     {
2003         GMX_RELEASE_ASSERT(link->a || link->index[cg_gl+1]+1 > link->nalloc_a,
2004                            "Inconsistent allocation of link");
2005         /* Add this charge group link */
2006         if (link->index[cg_gl+1]+1 > link->nalloc_a)
2007         {
2008             link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
2009             srenew(link->a, link->nalloc_a);
2010         }
2011         link->a[link->index[cg_gl+1]] = cg_gl_j;
2012         link->index[cg_gl+1]++;
2013     }
2014 }
2015
2016 static int *make_at2cg(t_block *cgs)
2017 {
2018     int *at2cg, cg, a;
2019
2020     snew(at2cg, cgs->index[cgs->nr]);
2021     for (cg = 0; cg < cgs->nr; cg++)
2022     {
2023         for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
2024         {
2025             at2cg[a] = cg;
2026         }
2027     }
2028
2029     return at2cg;
2030 }
2031
2032 t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd,
2033                                   cginfo_mb_t *cginfo_mb)
2034 {
2035     gmx_reverse_top_t  *rt;
2036     int                 mb, cg_offset, cg, cg_gl, a, aj, i, j, ftype, nral, nlink_mol, mol, ncgi;
2037     gmx_molblock_t     *molb;
2038     gmx_moltype_t      *molt;
2039     t_block            *cgs;
2040     t_blocka           *excls;
2041     int                *a2c;
2042     gmx_reverse_ilist_t ril;
2043     t_blocka           *link;
2044     cginfo_mb_t        *cgi_mb;
2045
2046     /* For each charge group make a list of other charge groups
2047      * in the system that a linked to it via bonded interactions
2048      * which are also stored in reverse_top.
2049      */
2050
2051     rt = dd->reverse_top;
2052
2053     snew(link, 1);
2054     snew(link->index, ncg_mtop(mtop)+1);
2055     link->nalloc_a = 0;
2056     link->a        = NULL;
2057
2058     link->index[0] = 0;
2059     cg_offset      = 0;
2060     ncgi           = 0;
2061     for (mb = 0; mb < mtop->nmolblock; mb++)
2062     {
2063         molb = &mtop->molblock[mb];
2064         if (molb->nmol == 0)
2065         {
2066             continue;
2067         }
2068         molt  = &mtop->moltype[molb->type];
2069         cgs   = &molt->cgs;
2070         excls = &molt->excls;
2071         a2c   = make_at2cg(cgs);
2072         /* Make a reverse ilist in which the interactions are linked
2073          * to all atoms, not only the first atom as in gmx_reverse_top.
2074          * The constraints are discarded here.
2075          */
2076         make_reverse_ilist(molt, NULL, FALSE, FALSE, FALSE, TRUE, &ril);
2077
2078         cgi_mb = &cginfo_mb[mb];
2079
2080         for (cg = 0; cg < cgs->nr; cg++)
2081         {
2082             cg_gl                = cg_offset + cg;
2083             link->index[cg_gl+1] = link->index[cg_gl];
2084             for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
2085             {
2086                 i = ril.index[a];
2087                 while (i < ril.index[a+1])
2088                 {
2089                     ftype = ril.il[i++];
2090                     nral  = NRAL(ftype);
2091                     /* Skip the ifunc index */
2092                     i++;
2093                     for (j = 0; j < nral; j++)
2094                     {
2095                         aj = ril.il[i+j];
2096                         if (a2c[aj] != cg)
2097                         {
2098                             check_link(link, cg_gl, cg_offset+a2c[aj]);
2099                         }
2100                     }
2101                     i += nral_rt(ftype);
2102                 }
2103                 if (rt->bExclRequired)
2104                 {
2105                     /* Exclusions always go both ways */
2106                     for (j = excls->index[a]; j < excls->index[a+1]; j++)
2107                     {
2108                         aj = excls->a[j];
2109                         if (a2c[aj] != cg)
2110                         {
2111                             check_link(link, cg_gl, cg_offset+a2c[aj]);
2112                         }
2113                     }
2114                 }
2115             }
2116             if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
2117             {
2118                 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
2119                 ncgi++;
2120             }
2121         }
2122         nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
2123
2124         cg_offset += cgs->nr;
2125
2126         destroy_reverse_ilist(&ril);
2127         sfree(a2c);
2128
2129         if (debug)
2130         {
2131             fprintf(debug, "molecule type '%s' %d cgs has %d cg links through bonded interac.\n", *molt->name, cgs->nr, nlink_mol);
2132         }
2133
2134         if (molb->nmol > 1)
2135         {
2136             /* Copy the data for the rest of the molecules in this block */
2137             link->nalloc_a += (molb->nmol - 1)*nlink_mol;
2138             srenew(link->a, link->nalloc_a);
2139             for (mol = 1; mol < molb->nmol; mol++)
2140             {
2141                 for (cg = 0; cg < cgs->nr; cg++)
2142                 {
2143                     cg_gl                = cg_offset + cg;
2144                     link->index[cg_gl+1] =
2145                         link->index[cg_gl+1-cgs->nr] + nlink_mol;
2146                     for (j = link->index[cg_gl]; j < link->index[cg_gl+1]; j++)
2147                     {
2148                         link->a[j] = link->a[j-nlink_mol] + cgs->nr;
2149                     }
2150                     if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
2151                         cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
2152                     {
2153                         SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
2154                         ncgi++;
2155                     }
2156                 }
2157                 cg_offset += cgs->nr;
2158             }
2159         }
2160     }
2161
2162     if (debug)
2163     {
2164         fprintf(debug, "Of the %d charge groups %d are linked via bonded interactions\n", ncg_mtop(mtop), ncgi);
2165     }
2166
2167     return link;
2168 }
2169
2170 static void bonded_cg_distance_mol(gmx_moltype_t *molt, int *at2cg,
2171                                    gmx_bool bBCheck, gmx_bool bExcl, rvec *cg_cm,
2172                                    real *r_2b, int *ft2b, int *a2_1, int *a2_2,
2173                                    real *r_mb, int *ftmb, int *am_1, int *am_2)
2174 {
2175     int       ftype, nral, i, j, ai, aj, cgi, cgj;
2176     t_ilist  *il;
2177     t_blocka *excls;
2178     real      r2_2b, r2_mb, rij2;
2179
2180     r2_2b = 0;
2181     r2_mb = 0;
2182     for (ftype = 0; ftype < F_NRE; ftype++)
2183     {
2184         if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE))
2185         {
2186             il   = &molt->ilist[ftype];
2187             nral = NRAL(ftype);
2188             if (nral > 1)
2189             {
2190                 for (i = 0; i < il->nr; i += 1+nral)
2191                 {
2192                     for (ai = 0; ai < nral; ai++)
2193                     {
2194                         cgi = at2cg[il->iatoms[i+1+ai]];
2195                         for (aj = 0; aj < nral; aj++)
2196                         {
2197                             cgj = at2cg[il->iatoms[i+1+aj]];
2198                             if (cgi != cgj)
2199                             {
2200                                 rij2 = distance2(cg_cm[cgi], cg_cm[cgj]);
2201                                 if (nral == 2 && rij2 > r2_2b)
2202                                 {
2203                                     r2_2b = rij2;
2204                                     *ft2b = ftype;
2205                                     *a2_1 = il->iatoms[i+1+ai];
2206                                     *a2_2 = il->iatoms[i+1+aj];
2207                                 }
2208                                 if (nral >  2 && rij2 > r2_mb)
2209                                 {
2210                                     r2_mb = rij2;
2211                                     *ftmb = ftype;
2212                                     *am_1 = il->iatoms[i+1+ai];
2213                                     *am_2 = il->iatoms[i+1+aj];
2214                                 }
2215                             }
2216                         }
2217                     }
2218                 }
2219             }
2220         }
2221     }
2222     if (bExcl)
2223     {
2224         excls = &molt->excls;
2225         for (ai = 0; ai < excls->nr; ai++)
2226         {
2227             cgi = at2cg[ai];
2228             for (j = excls->index[ai]; j < excls->index[ai+1]; j++)
2229             {
2230                 cgj = at2cg[excls->a[j]];
2231                 if (cgi != cgj)
2232                 {
2233                     rij2 = distance2(cg_cm[cgi], cg_cm[cgj]);
2234                     if (rij2 > r2_2b)
2235                     {
2236                         r2_2b = rij2;
2237                     }
2238                 }
2239             }
2240         }
2241     }
2242
2243     *r_2b = sqrt(r2_2b);
2244     *r_mb = sqrt(r2_mb);
2245 }
2246
2247 static void get_cgcm_mol(gmx_moltype_t *molt, gmx_ffparams_t *ffparams,
2248                          int ePBC, t_graph *graph, matrix box,
2249                          gmx_vsite_t *vsite,
2250                          rvec *x, rvec *xs, rvec *cg_cm)
2251 {
2252     int n, i;
2253
2254     if (ePBC != epbcNONE)
2255     {
2256         mk_mshift(NULL, graph, ePBC, box, x);
2257
2258         shift_x(graph, box, x, xs);
2259         /* By doing an extra mk_mshift the molecules that are broken
2260          * because they were e.g. imported from another software
2261          * will be made whole again. Such are the healing powers
2262          * of GROMACS.
2263          */
2264         mk_mshift(NULL, graph, ePBC, box, xs);
2265     }
2266     else
2267     {
2268         /* We copy the coordinates so the original coordinates remain
2269          * unchanged, just to be 100% sure that we do not affect
2270          * binary reproducibility of simulations.
2271          */
2272         n = molt->cgs.index[molt->cgs.nr];
2273         for (i = 0; i < n; i++)
2274         {
2275             copy_rvec(x[i], xs[i]);
2276         }
2277     }
2278
2279     if (vsite)
2280     {
2281         construct_vsites(vsite, xs, 0.0, NULL,
2282                          ffparams->iparams, molt->ilist,
2283                          epbcNONE, TRUE, NULL, NULL);
2284     }
2285
2286     calc_cgcm(NULL, 0, molt->cgs.nr, &molt->cgs, xs, cg_cm);
2287 }
2288
2289 static int have_vsite_molt(gmx_moltype_t *molt)
2290 {
2291     int      i;
2292     gmx_bool bVSite;
2293
2294     bVSite = FALSE;
2295     for (i = 0; i < F_NRE; i++)
2296     {
2297         if ((interaction_function[i].flags & IF_VSITE) &&
2298             molt->ilist[i].nr > 0)
2299         {
2300             bVSite = TRUE;
2301         }
2302     }
2303
2304     return bVSite;
2305 }
2306
2307 void dd_bonded_cg_distance(FILE *fplog,
2308                            gmx_mtop_t *mtop,
2309                            t_inputrec *ir, rvec *x, matrix box,
2310                            gmx_bool bBCheck,
2311                            real *r_2b, real *r_mb)
2312 {
2313     gmx_bool        bExclRequired;
2314     int             mb, at_offset, *at2cg, mol;
2315     t_graph         graph;
2316     gmx_vsite_t    *vsite;
2317     gmx_molblock_t *molb;
2318     gmx_moltype_t  *molt;
2319     rvec           *xs, *cg_cm;
2320     real            rmol_2b, rmol_mb;
2321     int             ft2b  = -1, a_2b_1 = -1, a_2b_2 = -1, ftmb = -1, a_mb_1 = -1, a_mb_2 = -1;
2322     int             ftm2b = -1, amol_2b_1 = -1, amol_2b_2 = -1, ftmmb = -1, amol_mb_1 = -1, amol_mb_2 = -1;
2323
2324     bExclRequired = IR_EXCL_FORCES(*ir);
2325
2326     vsite = init_vsite(mtop, NULL, TRUE);
2327
2328     *r_2b     = 0;
2329     *r_mb     = 0;
2330     at_offset = 0;
2331     for (mb = 0; mb < mtop->nmolblock; mb++)
2332     {
2333         molb = &mtop->molblock[mb];
2334         molt = &mtop->moltype[molb->type];
2335         if (molt->cgs.nr == 1 || molb->nmol == 0)
2336         {
2337             at_offset += molb->nmol*molt->atoms.nr;
2338         }
2339         else
2340         {
2341             if (ir->ePBC != epbcNONE)
2342             {
2343                 mk_graph_ilist(NULL, molt->ilist, 0, molt->atoms.nr, FALSE, FALSE,
2344                                &graph);
2345             }
2346
2347             at2cg = make_at2cg(&molt->cgs);
2348             snew(xs, molt->atoms.nr);
2349             snew(cg_cm, molt->cgs.nr);
2350             for (mol = 0; mol < molb->nmol; mol++)
2351             {
2352                 get_cgcm_mol(molt, &mtop->ffparams, ir->ePBC, &graph, box,
2353                              have_vsite_molt(molt) ? vsite : NULL,
2354                              x+at_offset, xs, cg_cm);
2355
2356                 bonded_cg_distance_mol(molt, at2cg, bBCheck, bExclRequired, cg_cm,
2357                                        &rmol_2b, &ftm2b, &amol_2b_1, &amol_2b_2,
2358                                        &rmol_mb, &ftmmb, &amol_mb_1, &amol_mb_2);
2359                 if (rmol_2b > *r_2b)
2360                 {
2361                     *r_2b  = rmol_2b;
2362                     ft2b   = ftm2b;
2363                     a_2b_1 = at_offset + amol_2b_1;
2364                     a_2b_2 = at_offset + amol_2b_2;
2365                 }
2366                 if (rmol_mb > *r_mb)
2367                 {
2368                     *r_mb  = rmol_mb;
2369                     ftmb   = ftmmb;
2370                     a_mb_1 = at_offset + amol_mb_1;
2371                     a_mb_2 = at_offset + amol_mb_2;
2372                 }
2373
2374                 at_offset += molt->atoms.nr;
2375             }
2376             sfree(cg_cm);
2377             sfree(xs);
2378             sfree(at2cg);
2379             if (ir->ePBC != epbcNONE)
2380             {
2381                 done_graph(&graph);
2382             }
2383         }
2384     }
2385
2386     /* We should have a vsite free routine, but here we can simply free */
2387     sfree(vsite);
2388
2389     if (fplog && (ft2b >= 0 || ftmb >= 0))
2390     {
2391         fprintf(fplog,
2392                 "Initial maximum inter charge-group distances:\n");
2393         if (ft2b >= 0)
2394         {
2395             fprintf(fplog,
2396                     "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2397                     *r_2b, interaction_function[ft2b].longname,
2398                     a_2b_1+1, a_2b_2+1);
2399         }
2400         if (ftmb >= 0)
2401         {
2402             fprintf(fplog,
2403                     "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2404                     *r_mb, interaction_function[ftmb].longname,
2405                     a_mb_1+1, a_mb_2+1);
2406         }
2407     }
2408 }