src/mdlib/domdec_top.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * This file is part of Gromacs        Copyright (c) 1991-2008
   5  * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation; either version 2
  10  * of the License, or (at your option) any later version.
  11  *
  12  * To help us fund GROMACS development, we humbly ask that you cite
  13  * the research papers on the package. Check out http://www.gromacs.org
  14  *
  15  * And Hey:
  16  * Gnomes, ROck Monsters And Chili Sauce
  17  */
  18
  19 #ifdef HAVE_CONFIG_H
  20 #include <config.h>
  21 #endif
  22
  23 #include <string.h>
  24 #include "typedefs.h"
  25 #include "smalloc.h"
  26 #include "domdec.h"
  27 #include "domdec_network.h"
  28 #include "names.h"
  29 #include "network.h"
  30 #include "vec.h"
  31 #include "pbc.h"
  32 #include "gmx_random.h"
  33 #include "topsort.h"
  34 #include "mtop_util.h"
  35 #include "mshift.h"
  36 #include "vsite.h"
  37 #include "gmx_ga2la.h"
  38
  39 typedef struct {
  40     int  *index;  /* Index for each atom into il                  */
  41     int  *il;     /* ftype|type|a0|...|an|ftype|...               */
  42 } gmx_reverse_ilist_t;
  43
  44 typedef struct {
  45     int  a_start;
  46     int  a_end;
  47     int  natoms_mol;
  48     int  type;
  49 } gmx_molblock_ind_t;
  50
  51 typedef struct gmx_reverse_top {
  52     bool bExclRequired; /* Do we require all exclusions to be assigned? */
  53     bool bConstr;       /* Are there constraints in this revserse top?  */
  54     bool bBCheck;       /* All bonded interactions have to be assigned? */
  55     bool bMultiCGmols;  /* Are the multi charge-group molecules?        */
  56     gmx_reverse_ilist_t *ril_mt; /* Reverse ilist for all moltypes      */
  57     int  ril_mt_tot_size;
  58     int  ilsort;        /* The sorting state of bondeds for free energy */
  59     gmx_molblock_ind_t *mbi;
  60
  61     /* Pointers only used for an error message */
  62     gmx_mtop_t     *err_top_global;
  63     gmx_localtop_t *err_top_local;
  64 } gmx_reverse_top_t;
  65
  66 static int nral_rt(int ftype)
  67 {
  68     /* Returns the number of atom entries for il in gmx_reverse_top_t */
  69     int nral;
  70
  71     nral = NRAL(ftype);
  72     if (interaction_function[ftype].flags & IF_VSITE)
  73     {
  74         /* With vsites the reverse topology contains
  75          * two extra entries for PBC.
  76          */
  77         nral += 2;
  78     }
  79
  80     return nral;
  81 }
  82
  83 static bool dd_check_ftype(int ftype,bool bBCheck,bool bConstr)
  84 {
  85     return (((interaction_function[ftype].flags & IF_BOND) &&
  86              !(interaction_function[ftype].flags & IF_VSITE) &&
  87              (bBCheck || !(interaction_function[ftype].flags & IF_LIMZERO))) ||
  88             ftype == F_SETTLE ||
  89             (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)));
  90 }
  91
  92 static void print_error_header(FILE *fplog,char *moltypename,int nprint)
  93 {
  94     fprintf(fplog, "\nMolecule type '%s'\n",moltypename);
  95     fprintf(stderr,"\nMolecule type '%s'\n",moltypename);
  96     fprintf(fplog,
  97             "the first %d missing interactions, except for exclusions:\n",
  98             nprint);
  99     fprintf(stderr,
 100             "the first %d missing interactions, except for exclusions:\n",
 101             nprint);
 102 }
 103
 104 static void print_missing_interactions_mb(FILE *fplog,t_commrec *cr,
 105                                           gmx_reverse_top_t *rt,
 106                                           char *moltypename,
 107                                           gmx_reverse_ilist_t *ril,
 108                                           int a_start,int a_end,
 109                                           int nat_mol,int nmol,
 110                                           t_idef *idef)
 111 {
 112     int nril_mol,*assigned,*gatindex;
 113     int ftype,ftype_j,nral,i,j_mol,j,k,a0,a0_mol,mol,a,a_gl;
 114     int nprint;
 115     t_ilist *il;
 116     t_iatom *ia;
 117     bool bFound;
 118
 119     nril_mol = ril->index[nat_mol];
 120     snew(assigned,nmol*nril_mol);
 121
 122     gatindex = cr->dd->gatindex;
 123     for(ftype=0; ftype<F_NRE; ftype++)
 124     {
 125         if (dd_check_ftype(ftype,rt->bBCheck,rt->bConstr))
 126         {
 127             nral = NRAL(ftype);
 128             il = &idef->il[ftype];
 129             ia = il->iatoms;
 130             for(i=0; i<il->nr; i+=1+nral)
 131             {
 132                 a0     = gatindex[ia[1]];
 133                 /* Check if this interaction is in
 134                  * the currently checked molblock.
 135                  */
 136                 if (a0 >= a_start && a0 < a_end)
 137                 {
 138                     mol    = (a0 - a_start)/nat_mol;
 139                     a0_mol = (a0 - a_start) - mol*nat_mol;
 140                     j_mol  = ril->index[a0_mol];
 141                     bFound = FALSE;
 142                     while (j_mol < ril->index[a0_mol+1] && !bFound)
 143                     {
 144                         j = mol*nril_mol + j_mol;
 145                         ftype_j = ril->il[j_mol];
 146                         /* Here we need to check if this interaction has
 147                          * not already been assigned, since we could have
 148                          * multiply defined interactions.
 149                          */
 150                         if (ftype == ftype_j && ia[0] == ril->il[j_mol+1] &&
 151                             assigned[j] == 0)
 152                         {
 153                             /* Check the atoms */
 154                             bFound = TRUE;
 155                             for(a=0; a<nral; a++)
 156                             {
 157                                 if (gatindex[ia[1+a]] !=
 158                                     a_start + mol*nat_mol + ril->il[j_mol+2+a])
 159                                 {
 160                                     bFound = FALSE;
 161                                 }
 162                             }
 163                             if (bFound)
 164                             {
 165                                 assigned[j] = 1;
 166                             }
 167                         }
 168                         j_mol += 2 + nral_rt(ftype_j);
 169                     }
 170                     if (!bFound)
 171                     {
 172                         gmx_incons("Some interactions seem to be assigned multiple times");
 173                     }
 174                 }
 175                 ia += 1 + nral;
 176             }
 177         }
 178     }
 179
 180     gmx_sumi(nmol*nril_mol,assigned,cr);
 181
 182     nprint = 10;
 183     i = 0;
 184     for(mol=0; mol<nmol; mol++)
 185     {
 186         j_mol = 0;
 187         while (j_mol < nril_mol)
 188         {
 189             ftype = ril->il[j_mol];
 190             nral  = NRAL(ftype);
 191             j = mol*nril_mol + j_mol;
 192             if (assigned[j] == 0 &&
 193                 !(interaction_function[ftype].flags & IF_VSITE))
 194             {
 195                 if (DDMASTER(cr->dd))
 196                 {
 197                     if (i == 0)
 198                     {
 199                         print_error_header(fplog,moltypename,nprint);
 200                     }
 201                     fprintf(fplog, "%20s atoms",
 202                             interaction_function[ftype].longname);
 203                     fprintf(stderr,"%20s atoms",
 204                             interaction_function[ftype].longname);
 205                     for(a=0; a<nral; a++) {
 206                         fprintf(fplog, "%5d",ril->il[j_mol+2+a]+1);
 207                         fprintf(stderr,"%5d",ril->il[j_mol+2+a]+1);
 208                     }
 209                     while (a < 4)
 210                     {
 211                         fprintf(fplog, "     ");
 212                         fprintf(stderr,"     ");
 213                         a++;
 214                     }
 215                     fprintf(fplog, " global");
 216                     fprintf(stderr," global");
 217                     for(a=0; a<nral; a++)
 218                     {
 219                         fprintf(fplog, "%6d",
 220                                 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
 221                         fprintf(stderr,"%6d",
 222                                 a_start+mol*nat_mol+ril->il[j_mol+2+a]+1);
 223                     }
 224                     fprintf(fplog, "\n");
 225                     fprintf(stderr,"\n");
 226                 }
 227                 i++;
 228                 if (i >= nprint)
 229                 {
 230                     break;
 231                 }
 232             }
 233             j_mol += 2 + nral_rt(ftype);
 234         }
 235     }
 236
 237     sfree(assigned);
 238 }
 239
 240 static void print_missing_interactions_atoms(FILE *fplog,t_commrec *cr,
 241                                              gmx_mtop_t *mtop,t_idef *idef)
 242 {
 243     int mb,a_start,a_end;
 244     gmx_molblock_t *molb;
 245     gmx_reverse_top_t *rt;
 246
 247     rt = cr->dd->reverse_top;
 248
 249     /* Print the atoms in the missing interactions per molblock */
 250     a_end = 0;
 251     for(mb=0; mb<mtop->nmolblock; mb++)
 252     {
 253         molb = &mtop->molblock[mb];
 254         a_start = a_end;
 255         a_end   = a_start + molb->nmol*molb->natoms_mol;
 256
 257         print_missing_interactions_mb(fplog,cr,rt,
 258                                       *(mtop->moltype[molb->type].name),
 259                                       &rt->ril_mt[molb->type],
 260                                       a_start,a_end,molb->natoms_mol,
 261                                       molb->nmol,
 262                                       idef);
 263     }
 264 }
 265
 266 void dd_print_missing_interactions(FILE *fplog,t_commrec *cr,int local_count,  gmx_mtop_t *top_global, t_state *state_local)
 267 {
 268     int  ndiff_tot,cl[F_NRE],n,ndiff,rest_global,rest_local;
 269     int  ftype,nral;
 270     char buf[STRLEN];
 271     gmx_domdec_t *dd;
 272     gmx_mtop_t     *err_top_global;
 273     gmx_localtop_t *err_top_local;
 274
 275     dd = cr->dd;
 276
 277     err_top_global = dd->reverse_top->err_top_global;
 278     err_top_local  = dd->reverse_top->err_top_local;
 279
 280     if (fplog)
 281     {
 282         fprintf(fplog,"\nNot all bonded interactions have been properly assigned to the domain decomposition cells\n");
 283         fflush(fplog);
 284     }
 285
 286     ndiff_tot = local_count - dd->nbonded_global;
 287
 288     for(ftype=0; ftype<F_NRE; ftype++)
 289     {
 290         nral = NRAL(ftype);
 291         cl[ftype] = err_top_local->idef.il[ftype].nr/(1+nral);
 292     }
 293
 294     gmx_sumi(F_NRE,cl,cr);
 295
 296     if (DDMASTER(dd))
 297     {
 298         fprintf(fplog,"\nA list of missing interactions:\n");
 299         fprintf(stderr,"\nA list of missing interactions:\n");
 300         rest_global = dd->nbonded_global;
 301         rest_local  = local_count;
 302         for(ftype=0; ftype<F_NRE; ftype++)
 303         {
 304             /* In the reverse and local top all constraints are merged
 305              * into F_CONSTR. So in the if statement we skip F_CONSTRNC
 306              * and add these constraints when doing F_CONSTR.
 307              */
 308             if (((interaction_function[ftype].flags & IF_BOND) &&
 309                  (dd->reverse_top->bBCheck
 310                   || !(interaction_function[ftype].flags & IF_LIMZERO)))
 311                 || ftype == F_SETTLE
 312                 || (dd->reverse_top->bConstr && ftype == F_CONSTR))
 313             {
 314                 nral = NRAL(ftype);
 315                 n = gmx_mtop_ftype_count(err_top_global,ftype);
 316                 if (ftype == F_CONSTR)
 317                 {
 318                     n += gmx_mtop_ftype_count(err_top_global,F_CONSTRNC);
 319                 }
 320                 ndiff = cl[ftype] - n;
 321                 if (ndiff != 0)
 322                 {
 323                     sprintf(buf,"%20s of %6d missing %6d",
 324                             interaction_function[ftype].longname,n,-ndiff);
 325                     fprintf(fplog,"%s\n",buf);
 326                     fprintf(stderr,"%s\n",buf);
 327                 }
 328                 rest_global -= n;
 329                 rest_local  -= cl[ftype];
 330             }
 331         }
 332
 333         ndiff = rest_local - rest_global;
 334         if (ndiff != 0)
 335         {
 336             sprintf(buf,"%20s of %6d missing %6d","exclusions",
 337                     rest_global,-ndiff);
 338             fprintf(fplog,"%s\n",buf);
 339             fprintf(stderr,"%s\n",buf);
 340         }
 341     }
 342
 343     print_missing_interactions_atoms(fplog,cr,err_top_global,
 344                                      &err_top_local->idef);
 345     write_dd_pdb("dd_dump_err",0,"dump",top_global,cr,
 346                  -1,state_local->x,state_local->box);
 347     if (DDMASTER(dd))
 348     {
 349         if (ndiff_tot > 0)
 350         {
 351             gmx_incons("One or more interactions were multiple assigned in the domain decompostion");
 352         }
 353         else
 354         {
 355             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));
 356         }
 357     }
 358 }
 359
 360 static void global_atomnr_to_moltype_ind(gmx_molblock_ind_t *mbi,int i_gl,
 361                                          int *mb,int *mt,int *mol,int *i_mol)
 362 {
 363   int molb;
 364
 365   *mb = 0;
 366   while (i_gl >= mbi->a_end) {
 367     (*mb)++;
 368     mbi++;
 369   }
 370
 371   *mt    = mbi->type;
 372   *mol   = (i_gl - mbi->a_start) / mbi->natoms_mol;
 373   *i_mol = (i_gl - mbi->a_start) - (*mol)*mbi->natoms_mol;
 374 }
 375
 376 static int count_excls(t_block *cgs,t_blocka *excls,int *n_intercg_excl)
 377 {
 378     int n,n_inter,cg,at0,at1,at,excl,atj;
 379
 380     n = 0;
 381     *n_intercg_excl = 0;
 382     for(cg=0; cg<cgs->nr; cg++)
 383     {
 384         at0 = cgs->index[cg];
 385         at1 = cgs->index[cg+1];
 386         for(at=at0; at<at1; at++)
 387         {
 388             for(excl=excls->index[at]; excl<excls->index[at+1]; excl++)
 389             {
 390                 atj = excls->a[excl];
 391                 if (atj > at)
 392                 {
 393                     n++;
 394                     if (atj < at0 || atj >= at1)
 395                     {
 396                         (*n_intercg_excl)++;
 397                     }
 398                 }
 399             }
 400         }
 401     }
 402
 403     return n;
 404 }
 405
 406 static int low_make_reverse_ilist(t_ilist *il_mt,t_atom *atom,
 407                                   int **vsite_pbc,
 408                                   int *count,
 409                                   bool bConstr,bool bBCheck,
 410                                   int *r_index,int *r_il,
 411                                   bool bLinkToAllAtoms,
 412                                   bool bAssign)
 413 {
 414     int  ftype,nral,i,j,nlink,link;
 415     t_ilist *il;
 416     t_iatom *ia;
 417     atom_id a;
 418     int  nint;
 419     bool bVSite;
 420
 421     nint = 0;
 422     for(ftype=0; ftype<F_NRE; ftype++)
 423     {
 424         if ((interaction_function[ftype].flags & (IF_BOND | IF_VSITE)) ||
 425             ftype == F_SETTLE ||
 426             (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC))) {
 427             bVSite = (interaction_function[ftype].flags & IF_VSITE);
 428             nral = NRAL(ftype);
 429             il = &il_mt[ftype];
 430             ia  = il->iatoms;
 431             for(i=0; i<il->nr; i+=1+nral)
 432             {
 433                 ia = il->iatoms + i;
 434                 if (bLinkToAllAtoms)
 435                 {
 436                     if (bVSite)
 437                     {
 438                         /* We don't need the virtual sites for the cg-links */
 439                         nlink = 0;
 440                     }
 441                     else
 442                     {
 443                         nlink = nral;
 444                     }
 445                 }
 446                 else
 447                 {
 448                     /* Couple to the first atom in the interaction */
 449                     nlink = 1;
 450                 }
 451                 for(link=0; link<nlink; link++)
 452                 {
 453                     a = ia[1+link];
 454                     if (bAssign)
 455                     {
 456                         r_il[r_index[a]+count[a]] =
 457                             (ftype == F_CONSTRNC ? F_CONSTR : ftype);
 458                         r_il[r_index[a]+count[a]+1] = ia[0];
 459                         for(j=1; j<1+nral; j++)
 460                         {
 461                             /* Store the molecular atom number */
 462                             r_il[r_index[a]+count[a]+1+j] = ia[j];
 463                         }
 464                     }
 465                     if (interaction_function[ftype].flags & IF_VSITE)
 466                     {
 467                         if (bAssign)
 468                         {
 469                             /* Add an entry to iatoms for storing
 470                              * which of the constructing atoms are
 471                              * vsites again.
 472                              */
 473                             r_il[r_index[a]+count[a]+2+nral] = 0;
 474                             for(j=2; j<1+nral; j++)
 475                             {
 476                                 if (atom[ia[j]].ptype == eptVSite)
 477                                 {
 478                                     r_il[r_index[a]+count[a]+2+nral] |= (2<<j);
 479                                 }
 480                             }
 481                             /* Store vsite pbc atom in a second extra entry */
 482                             r_il[r_index[a]+count[a]+2+nral+1] =
 483                                 (vsite_pbc ? vsite_pbc[ftype-F_VSITE2][i/(1+nral)] : -2);
 484                         }
 485                     }
 486                     else
 487                     {
 488                         /* We do not count vsites since they are always
 489                          * uniquely assigned and can be assigned
 490                          * to multiple nodes with recursive vsites.
 491                          */
 492                         if (bBCheck ||
 493                             !(interaction_function[ftype].flags & IF_LIMZERO))
 494                         {
 495                             nint++;
 496                         }
 497                     }
 498                     count[a] += 2 + nral_rt(ftype);
 499                 }
 500             }
 501         }
 502     }
 503
 504     return nint;
 505 }
 506
 507 static int make_reverse_ilist(gmx_moltype_t *molt,
 508                               int **vsite_pbc,
 509                               bool bConstr,bool bBCheck,
 510                               bool bLinkToAllAtoms,
 511                               gmx_reverse_ilist_t *ril_mt)
 512 {
 513     int nat_mt,*count,i,nint_mt;
 514
 515     /* Count the interactions */
 516     nat_mt = molt->atoms.nr;
 517     snew(count,nat_mt);
 518     low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
 519                            count,
 520                            bConstr,bBCheck,NULL,NULL,
 521                            bLinkToAllAtoms,FALSE);
 522
 523     snew(ril_mt->index,nat_mt+1);
 524     ril_mt->index[0] = 0;
 525     for(i=0; i<nat_mt; i++)
 526     {
 527         ril_mt->index[i+1] = ril_mt->index[i] + count[i];
 528         count[i] = 0;
 529     }
 530     snew(ril_mt->il,ril_mt->index[nat_mt]);
 531
 532     /* Store the interactions */
 533     nint_mt =
 534         low_make_reverse_ilist(molt->ilist,molt->atoms.atom,vsite_pbc,
 535                                count,
 536                                bConstr,bBCheck,
 537                                ril_mt->index,ril_mt->il,
 538                                bLinkToAllAtoms,TRUE);
 539
 540     sfree(count);
 541
 542     return nint_mt;
 543 }
 544
 545 static void destroy_reverse_ilist(gmx_reverse_ilist_t *ril)
 546 {
 547     sfree(ril->index);
 548     sfree(ril->il);
 549 }
 550
 551 static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop,bool bFE,
 552                                            int ***vsite_pbc_molt,
 553                                            bool bConstr,
 554                                            bool bBCheck,int *nint)
 555 {
 556     int mt,i,mb;
 557     gmx_reverse_top_t *rt;
 558     int *nint_mt;
 559     gmx_moltype_t *molt;
 560
 561     snew(rt,1);
 562
 563     /* Should we include constraints (for SHAKE) in rt? */
 564     rt->bConstr = bConstr;
 565     rt->bBCheck = bBCheck;
 566
 567     rt->bMultiCGmols = FALSE;
 568     snew(nint_mt,mtop->nmoltype);
 569     snew(rt->ril_mt,mtop->nmoltype);
 570     rt->ril_mt_tot_size = 0;
 571     for(mt=0; mt<mtop->nmoltype; mt++)
 572     {
 573         molt = &mtop->moltype[mt];
 574         if (molt->cgs.nr > 1)
 575         {
 576             rt->bMultiCGmols = TRUE;
 577         }
 578
 579         /* Make the atom to interaction list for this molecule type */
 580         nint_mt[mt] =
 581             make_reverse_ilist(molt,vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
 582                                rt->bConstr,rt->bBCheck,FALSE,
 583                                &rt->ril_mt[mt]);
 584
 585         rt->ril_mt_tot_size += rt->ril_mt[mt].index[molt->atoms.nr];
 586     }
 587     if (debug)
 588     {
 589         fprintf(debug,"The total size of the atom to interaction index is %d integers\n",rt->ril_mt_tot_size);
 590     }
 591
 592     *nint = 0;
 593     for(mb=0; mb<mtop->nmolblock; mb++)
 594     {
 595         *nint += mtop->molblock[mb].nmol*nint_mt[mtop->molblock[mb].type];
 596     }
 597     sfree(nint_mt);
 598
 599     if (bFE && gmx_mtop_bondeds_free_energy(mtop))
 600     {
 601         rt->ilsort = ilsortFE_UNSORTED;
 602     }
 603     else {
 604         rt->ilsort = ilsortNO_FE;
 605     }
 606
 607     /* Make a molblock index for fast searching */
 608     snew(rt->mbi,mtop->nmolblock);
 609     i = 0;
 610     for(mb=0; mb<mtop->nmolblock; mb++)
 611     {
 612         rt->mbi[mb].a_start    = i;
 613         i += mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
 614         rt->mbi[mb].a_end      = i;
 615         rt->mbi[mb].natoms_mol = mtop->molblock[mb].natoms_mol;
 616         rt->mbi[mb].type       = mtop->molblock[mb].type;
 617     }
 618
 619     return rt;
 620 }
 621
 622 void dd_make_reverse_top(FILE *fplog,
 623                          gmx_domdec_t *dd,gmx_mtop_t *mtop,
 624                          gmx_vsite_t *vsite,gmx_constr_t constr,
 625                          t_inputrec *ir,bool bBCheck)
 626 {
 627     int mb,natoms,n_recursive_vsite,nexcl,nexcl_icg,a;
 628     gmx_molblock_t *molb;
 629     gmx_moltype_t *molt;
 630
 631     if (fplog)
 632     {
 633         fprintf(fplog,"\nLinking all bonded interactions to atoms\n");
 634     }
 635
 636     dd->reverse_top = make_reverse_top(mtop,ir->efep!=efepNO,
 637                                        vsite ? vsite->vsite_pbc_molt : NULL,
 638                                        !dd->bInterCGcons,
 639                                        bBCheck,&dd->nbonded_global);
 640
 641     if (dd->reverse_top->ril_mt_tot_size >= 200000 &&
 642         mtop->mols.nr > 1 &&
 643         mtop->nmolblock == 1 && mtop->molblock[0].nmol == 1)
 644     {
 645         /* mtop comes from a pre Gromacs 4 tpr file */
 646         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";
 647         if (fplog)
 648         {
 649             fprintf(fplog,"\n%s\n\n",note);
 650         }
 651         if (DDMASTER(dd))
 652         {
 653             fprintf(stderr,"\n%s\n\n",note);
 654         }
 655     }
 656
 657     dd->reverse_top->bExclRequired = IR_EXCL_FORCES(*ir);
 658
 659     nexcl = 0;
 660     dd->n_intercg_excl = 0;
 661     for(mb=0; mb<mtop->nmolblock; mb++)
 662     {
 663         molb = &mtop->molblock[mb];
 664         molt = &mtop->moltype[molb->type];
 665         nexcl += molb->nmol*count_excls(&molt->cgs,&molt->excls,&nexcl_icg);
 666         dd->n_intercg_excl += molb->nmol*nexcl_icg;
 667     }
 668     if (dd->reverse_top->bExclRequired)
 669     {
 670         dd->nbonded_global += nexcl;
 671         if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0 && fplog)
 672         {
 673             fprintf(fplog,"There are %d inter charge-group exclusions,\n"
 674                     "will use an extra communication step for exclusion forces for %s\n",
 675                     dd->n_intercg_excl,eel_names[ir->coulombtype]);
 676         }
 677     }
 678
 679     natoms = mtop->natoms;
 680
 681     if (vsite && vsite->n_intercg_vsite > 0)
 682     {
 683         if (fplog)
 684         {
 685             fprintf(fplog,"There are %d inter charge-group virtual sites,\n"
 686                     "will an extra communication step for selected coordinates and forces\n",
 687               vsite->n_intercg_vsite);
 688         }
 689         init_domdec_vsites(dd,natoms);
 690     }
 691
 692     if (dd->bInterCGcons)
 693     {
 694         init_domdec_constraints(dd,natoms,mtop,constr);
 695     }
 696     if (fplog)
 697     {
 698         fprintf(fplog,"\n");
 699     }
 700 }
 701
 702 static inline void add_ifunc(int nral,t_iatom *tiatoms,t_ilist *il)
 703 {
 704     t_iatom *liatoms;
 705     int     k;
 706
 707     if (il->nr+1+nral > il->nalloc)
 708     {
 709         il->nalloc += over_alloc_large(il->nr+1+nral);
 710         srenew(il->iatoms,il->nalloc);
 711     }
 712     liatoms = il->iatoms + il->nr;
 713     for(k=0; k<=nral; k++)
 714     {
 715         liatoms[k] = tiatoms[k];
 716     }
 717     il->nr += 1 + nral;
 718 }
 719
 720 static void add_posres(int mol,int a_mol,gmx_molblock_t *molb,
 721                        t_iatom *iatoms,t_idef *idef)
 722 {
 723     int n,a_molb;
 724     t_iparams *ip;
 725
 726     /* This position restraint has not been added yet,
 727      * so it's index is the current number of position restraints.
 728      */
 729     n = idef->il[F_POSRES].nr/2;
 730     if (n >= idef->iparams_posres_nalloc)
 731     {
 732         idef->iparams_posres_nalloc = over_alloc_dd(n);
 733         srenew(idef->iparams_posres,idef->iparams_posres_nalloc);
 734     }
 735     ip = &idef->iparams_posres[n];
 736     /* Copy the force constants */
 737     *ip = idef->iparams[iatoms[0]];
 738
 739     /* Get the position restriant coordinats from the molblock */
 740     a_molb = mol*molb->natoms_mol + a_mol;
 741     if (a_molb >= molb->nposres_xA)
 742     {
 743         gmx_incons("Not enough position restraint coordinates");
 744     }
 745     ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
 746     ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
 747     ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
 748     if (molb->nposres_xB > 0)
 749     {
 750         ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
 751         ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
 752         ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
 753     }
 754     else
 755     {
 756         ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
 757         ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
 758         ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
 759     }
 760     /* Set the parameter index for idef->iparams_posre */
 761     iatoms[0] = n;
 762 }
 763
 764 static void add_vsite(gmx_ga2la_t ga2la,int *index,int *rtil,
 765                       int ftype,int nral,
 766                       bool bHomeA,int a,int a_gl,int a_mol,
 767                       t_iatom *iatoms,
 768                       t_idef *idef,int **vsite_pbc,int *vsite_pbc_nalloc)
 769 {
 770     int  k,ak_gl,vsi,pbc_a_mol;
 771     t_iatom tiatoms[1+MAXATOMLIST],*iatoms_r;
 772     int  j,ftype_r,nral_r;
 773
 774     /* Copy the type */
 775     tiatoms[0] = iatoms[0];
 776
 777     if (bHomeA)
 778     {
 779         /* We know the local index of the first atom */
 780         tiatoms[1] = a;
 781     }
 782     else
 783     {
 784         /* Convert later in make_local_vsites */
 785         tiatoms[1] = -a_gl - 1;
 786     }
 787
 788     for(k=2; k<1+nral; k++)
 789     {
 790         ak_gl = a_gl + iatoms[k] - a_mol;
 791         if (!ga2la_home(ga2la,ak_gl,&tiatoms[k]))
 792         {
 793             /* Copy the global index, convert later in make_local_vsites */
 794             tiatoms[k] = -(ak_gl + 1);
 795         }
 796     }
 797
 798     /* Add this interaction to the local topology */
 799     add_ifunc(nral,tiatoms,&idef->il[ftype]);
 800     if (vsite_pbc)
 801     {
 802         vsi = idef->il[ftype].nr/(1+nral) - 1;
 803         if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
 804         {
 805             vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
 806             srenew(vsite_pbc[ftype-F_VSITE2],vsite_pbc_nalloc[ftype-F_VSITE2]);
 807         }
 808         if (bHomeA)
 809         {
 810             pbc_a_mol = iatoms[1+nral+1];
 811             if (pbc_a_mol < 0)
 812             {
 813                 /* The pbc flag is one of the following two options:
 814                  * -2: vsite and all constructing atoms are within the same cg, no pbc
 815                  * -1: vsite and its first constructing atom are in the same cg, do pbc
 816                  */
 817                 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
 818             }
 819             else
 820             {
 821                 /* Set the pbc atom for this vsite so we can make its pbc
 822                  * identical to the rest of the atoms in its charge group.
 823                  * Since the order of the atoms does not change within a charge
 824                  * group, we do not need the global to local atom index.
 825                  */
 826                 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
 827             }
 828         }
 829         else
 830         {
 831             /* This vsite is non-home (required for recursion),
 832              * and therefore there is no charge group to match pbc with.
 833              * But we always turn on full_pbc to assure that higher order
 834              * recursion works correctly.
 835              */
 836             vsite_pbc[ftype-F_VSITE2][vsi] = -1;
 837         }
 838     }
 839
 840     if (iatoms[1+nral])
 841     {
 842         /* Check for recursion */
 843         for(k=2; k<1+nral; k++)
 844         {
 845             if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
 846             {
 847                 /* This construction atoms is a vsite and not a home atom */
 848                 if (gmx_debug_at)
 849                 {
 850                     fprintf(debug,"Constructing atom %d of vsite atom %d is a vsite and non-home\n",iatoms[k]+1,a_mol+1);
 851                 }
 852                 /* Find the vsite construction */
 853
 854                 /* Check all interactions assigned to this atom */
 855                 j = index[iatoms[k]];
 856                 while (j < index[iatoms[k]+1])
 857                 {
 858                     ftype_r = rtil[j++];
 859                     nral_r = NRAL(ftype_r);
 860                     if (interaction_function[ftype_r].flags & IF_VSITE)
 861                     {
 862                         /* Add this vsite (recursion) */
 863                         add_vsite(ga2la,index,rtil,ftype_r,nral_r,
 864                                   FALSE,-1,a_gl+iatoms[k]-iatoms[1],iatoms[k],
 865                                   rtil+j,idef,vsite_pbc,vsite_pbc_nalloc);
 866                         j += 1 + nral_r + 2;
 867                     }
 868                     else
 869                     {
 870                         j += 1 + nral_r;
 871                     }
 872                 }
 873             }
 874         }
 875     }
 876 }
 877
 878 static void make_la2lc(gmx_domdec_t *dd)
 879 {
 880     int *cgindex,*la2lc,cg,a;
 881
 882     cgindex = dd->cgindex;
 883
 884     if (dd->nat_tot > dd->la2lc_nalloc)
 885     {
 886         dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
 887         snew(dd->la2lc,dd->la2lc_nalloc);
 888     }
 889     la2lc = dd->la2lc;
 890
 891     /* Make the local atom to local cg index */
 892     for(cg=0; cg<dd->ncg_tot; cg++)
 893     {
 894         for(a=cgindex[cg]; a<cgindex[cg+1]; a++)
 895         {
 896             la2lc[a] = cg;
 897         }
 898     }
 899 }
 900
 901 static real dd_dist2(t_pbc *pbc_null,rvec *cg_cm,const int *la2lc,int i,int j)
 902 {
 903     rvec dx;
 904
 905     if (pbc_null)
 906     {
 907         pbc_dx_aiuc(pbc_null,cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
 908     }
 909     else
 910     {
 911         rvec_sub(cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
 912     }
 913
 914     return norm2(dx);
 915 }
 916
 917 static int make_local_bondeds(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
 918                               gmx_molblock_t *molb,
 919                               bool bRCheckMB,ivec rcheck,bool bRCheck2B,
 920                               real rc,
 921                               int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
 922                               t_idef *idef,gmx_vsite_t *vsite)
 923 {
 924     int nzone,nizone,ic,la0,la1,i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,d,k;
 925     int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
 926     t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
 927     bool bBCheck,bUse,bLocal;
 928     real rc2;
 929     ivec k_zero,k_plus;
 930     gmx_ga2la_t ga2la;
 931     int  a_loc;
 932     int  kc;
 933     gmx_domdec_ns_ranges_t *izone;
 934     gmx_reverse_top_t *rt;
 935     gmx_molblock_ind_t *mbi;
 936     int nbonded_local;
 937
 938     nzone  = zones->n;
 939     nizone = zones->nizone;
 940     izone  = zones->izone;
 941
 942     rc2 = rc*rc;
 943
 944     if (vsite && vsite->n_intercg_vsite > 0)
 945     {
 946         vsite_pbc        = vsite->vsite_pbc_loc;
 947         vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
 948     }
 949     else
 950     {
 951         vsite_pbc        = NULL;
 952         vsite_pbc_nalloc = NULL;
 953     }
 954
 955     rt = dd->reverse_top;
 956
 957     bBCheck = rt->bBCheck;
 958
 959     /* Clear the counts */
 960     for(ftype=0; ftype<F_NRE; ftype++)
 961     {
 962         idef->il[ftype].nr = 0;
 963     }
 964     nbonded_local = 0;
 965
 966     mbi = rt->mbi;
 967
 968     ga2la = dd->ga2la;
 969
 970     for(ic=0; ic<nzone; ic++)
 971     {
 972         la0 = dd->cgindex[zones->cg_range[ic]];
 973         la1 = dd->cgindex[zones->cg_range[ic+1]];
 974         for(i=la0; i<la1; i++)
 975         {
 976             /* Get the global atom number */
 977             i_gl = dd->gatindex[i];
 978             global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
 979             /* Check all interactions assigned to this atom */
 980             index = rt->ril_mt[mt].index;
 981             rtil  = rt->ril_mt[mt].il;
 982             j = index[i_mol];
 983             while (j < index[i_mol+1])
 984             {
 985                 ftype  = rtil[j++];
 986                 iatoms = rtil + j;
 987                 nral = NRAL(ftype);
 988                 if (interaction_function[ftype].flags & IF_VSITE)
 989                 {
 990                     /* The vsite construction goes where the vsite itself is */
 991                     if (ic == 0)
 992                     {
 993                         add_vsite(dd->ga2la,index,rtil,ftype,nral,
 994                                   TRUE,i,i_gl,i_mol,
 995                                   iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
 996                     }
 997                     j += 1 + nral + 2;
 998                 }
 999                 else
1000                 {
1001                     /* Copy the type */
1002                     tiatoms[0] = iatoms[0];
1003
1004                     if (nral == 1)
1005                     {
1006                         /* Assign single-body interactions to the home zone */
1007                         if (ic == 0)
1008                         {
1009                             bUse = TRUE;
1010                             tiatoms[1] = i;
1011                             if (ftype == F_POSRES)
1012                             {
1013                                 add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1014                             }
1015                         }
1016                         else
1017                         {
1018                             bUse = FALSE;
1019                         }
1020                     }
1021                     else if (nral == 2)
1022                     {
1023                         /* This is a two-body interaction, we can assign
1024                          * analogous to the non-bonded assignments.
1025                          */
1026                         if (!ga2la_get(ga2la,i_gl+iatoms[2]-i_mol,&a_loc,&kc))
1027                         {
1028                             bUse = FALSE;
1029                         }
1030                         else
1031                         {
1032                             if (kc >= nzone)
1033                             {
1034                                 kc -= nzone;
1035                             }
1036                             /* Check zone interaction assignments */
1037                             bUse = ((ic < nizone && ic <= kc &&
1038                                      izone[ic].j0 <= kc && kc < izone[ic].j1) ||
1039                                     (kc < nizone && ic >  kc &&
1040                                      izone[kc].j0 <= ic && ic < izone[kc].j1));
1041                             if (bUse)
1042                             {
1043                                 tiatoms[1] = i;
1044                                 tiatoms[2] = a_loc;
1045                                 /* If necessary check the cgcm distance */
1046                                 if (bRCheck2B &&
1047                                     dd_dist2(pbc_null,cg_cm,la2lc,
1048                                              tiatoms[1],tiatoms[2]) >= rc2)
1049                                 {
1050                                     bUse = FALSE;
1051                                 }
1052                             }
1053                         }
1054                     }
1055                     else
1056                     {
1057                         /* Assign this multi-body bonded interaction to
1058                          * the local node if we have all the atoms involved
1059                          * (local or communicated) and the minimum zone shift
1060                          * in each dimension is zero, for dimensions
1061                          * with 2 DD cells an extra check may be necessary.
1062                          */
1063                         bUse = TRUE;
1064                         clear_ivec(k_zero);
1065                         clear_ivec(k_plus);
1066                         for(k=1; k<=nral && bUse; k++)
1067                         {
1068                             bLocal = ga2la_get(ga2la,i_gl+iatoms[k]-i_mol,
1069                                                &a_loc,&kc);
1070                             if (!bLocal || kc >= zones->n)
1071                             {
1072                                 /* We do not have this atom of this interaction
1073                                  * locally, or it comes from more than one cell
1074                                  * away.
1075                                  */
1076                                 bUse = FALSE;
1077                             }
1078                             else
1079                             {
1080                                 tiatoms[k] = a_loc;
1081                                 for(d=0; d<DIM; d++)
1082                                 {
1083                                     if (zones->shift[kc][d] == 0)
1084                                     {
1085                                         k_zero[d] = k;
1086                                     }
1087                                     else
1088                                     {
1089                                         k_plus[d] = k;
1090                                     }
1091                                 }
1092                             }
1093                         }
1094                         bUse = (bUse &&
1095                                 k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1096                         if (bRCheckMB)
1097                         {
1098                             for(d=0; (d<DIM && bUse); d++)
1099                             {
1100                                 /* Check if the cg_cm distance falls within
1101                                  * the cut-off to avoid possible multiple
1102                                  * assignments of bonded interactions.
1103                                  */
1104                                 if (rcheck[d] &&
1105                                     k_plus[d] &&
1106                                     dd_dist2(pbc_null,cg_cm,la2lc,
1107                                              tiatoms[k_zero[d]],tiatoms[k_plus[d]]) >= rc2)
1108                                 {
1109                                     bUse = FALSE;
1110                                 }
1111                             }
1112                         }
1113                     }
1114                     if (bUse)
1115                     {
1116                         /* Add this interaction to the local topology */
1117                         add_ifunc(nral,tiatoms,&idef->il[ftype]);
1118                         /* Sum so we can check in global_stat
1119                          * if we have everything.
1120                          */
1121                         if (bBCheck ||
1122                             !(interaction_function[ftype].flags & IF_LIMZERO))
1123                         {
1124                             nbonded_local++;
1125                         }
1126                     }
1127                     j += 1 + nral;
1128                 }
1129             }
1130         }
1131     }
1132
1133     return nbonded_local;
1134 }
1135
1136 static int make_local_bondeds_intracg(gmx_domdec_t *dd,gmx_molblock_t *molb,
1137                                       t_idef *idef,gmx_vsite_t *vsite)
1138 {
1139     int i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,k;
1140     int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
1141     t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
1142     gmx_reverse_top_t *rt;
1143     gmx_molblock_ind_t *mbi;
1144     int nbonded_local;
1145
1146     if (vsite && vsite->n_intercg_vsite > 0)
1147     {
1148         vsite_pbc        = vsite->vsite_pbc_loc;
1149         vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1150     }
1151     else
1152     {
1153         vsite_pbc        = NULL;
1154         vsite_pbc_nalloc = NULL;
1155     }
1156
1157     /* Clear the counts */
1158     for(ftype=0; ftype<F_NRE; ftype++)
1159     {
1160         idef->il[ftype].nr = 0;
1161     }
1162     nbonded_local = 0;
1163
1164     rt = dd->reverse_top;
1165
1166     if (rt->ril_mt_tot_size == 0)
1167     {
1168         /* There are no interactions to assign */
1169         return nbonded_local;
1170     }
1171
1172     mbi = rt->mbi;
1173
1174     for(i=0; i<dd->nat_home; i++)
1175     {
1176         /* Get the global atom number */
1177         i_gl = dd->gatindex[i];
1178         global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
1179         /* Check all interactions assigned to this atom */
1180         index = rt->ril_mt[mt].index;
1181         rtil  = rt->ril_mt[mt].il;
1182         /* Check all interactions assigned to this atom */
1183         j = index[i_mol];
1184         while (j < index[i_mol+1])
1185         {
1186             ftype  = rtil[j++];
1187             iatoms = rtil + j;
1188             nral = NRAL(ftype);
1189             if (interaction_function[ftype].flags & IF_VSITE)
1190             {
1191                 /* The vsite construction goes where the vsite itself is */
1192                 add_vsite(dd->ga2la,index,rtil,ftype,nral,
1193                           TRUE,i,i_gl,i_mol,
1194                           iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1195                 j += 1 + nral + 2;
1196             }
1197             else
1198             {
1199                 /* Copy the type */
1200                 tiatoms[0] = iatoms[0];
1201                 tiatoms[1] = i;
1202                 for(k=2; k<=nral; k++)
1203                 {
1204                     tiatoms[k] = i + iatoms[k] - iatoms[1];
1205                 }
1206                 if (ftype == F_POSRES)
1207                 {
1208                     add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1209                 }
1210                 /* Add this interaction to the local topology */
1211                 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1212                 /* Sum so we can check in global_stat if we have everything */
1213                 nbonded_local++;
1214                 j += 1 + nral;
1215             }
1216         }
1217     }
1218
1219     return nbonded_local;
1220 }
1221
1222 static int make_local_exclusions(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1223                                  gmx_mtop_t *mtop,
1224                                  bool bRCheck,real rc,
1225                                  int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1226                                  t_forcerec *fr,
1227                                  t_blocka *lexcls)
1228 {
1229     int  nizone,n,count,ic,jla0,jla1,jla;
1230     int  cg,la0,la1,la,a_gl,mb,mt,mol,a_mol,j,aj_mol;
1231     t_blocka *excls;
1232     gmx_ga2la_t ga2la;
1233     int  a_loc;
1234     int  cell;
1235     gmx_molblock_ind_t *mbi;
1236     real rc2;
1237
1238     /* Since for RF and PME we need to loop over the exclusions
1239      * we should store each exclusion only once. This is done
1240      * using the same zone scheme as used for neighbor searching.
1241      * The exclusions involving non-home atoms are stored only
1242      * one way: atom j is in the excl list of i only for j > i,
1243      * where i and j are local atom numbers.
1244      */
1245
1246     lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1247     if (lexcls->nr+1 > lexcls->nalloc_index)
1248     {
1249         lexcls->nalloc_index = over_alloc_dd(lexcls->nr)+1;
1250         srenew(lexcls->index,lexcls->nalloc_index);
1251     }
1252
1253     mbi = dd->reverse_top->mbi;
1254
1255     ga2la = dd->ga2la;
1256
1257     rc2 = rc*rc;
1258
1259     if (dd->n_intercg_excl)
1260     {
1261         nizone = zones->nizone;
1262     }
1263     else
1264     {
1265         nizone = 1;
1266     }
1267     n = 0;
1268     count = 0;
1269     for(ic=0; ic<nizone; ic++)
1270     {
1271         jla0 = dd->cgindex[zones->izone[ic].jcg0];
1272         jla1 = dd->cgindex[zones->izone[ic].jcg1];
1273         for(cg=zones->cg_range[ic]; cg<zones->cg_range[ic+1]; cg++)
1274         {
1275             /* Here we assume the number of exclusions in one charge group
1276              * is never larger than 1000.
1277              */
1278             if (n+1000 > lexcls->nalloc_a)
1279             {
1280                 lexcls->nalloc_a = over_alloc_large(n+1000);
1281                 srenew(lexcls->a,lexcls->nalloc_a);
1282             }
1283             la0 = dd->cgindex[cg];
1284             la1 = dd->cgindex[cg+1];
1285             if (GET_CGINFO_EXCL_INTER(fr->cginfo[cg]) ||
1286                 !GET_CGINFO_EXCL_INTRA(fr->cginfo[cg]))
1287             {
1288                 /* Copy the exclusions from the global top */
1289                 for(la=la0; la<la1; la++) {
1290                     lexcls->index[la] = n;
1291                     a_gl = dd->gatindex[la];
1292                     global_atomnr_to_moltype_ind(mbi,a_gl,&mb,&mt,&mol,&a_mol);
1293                     excls = &mtop->moltype[mt].excls;
1294                     for(j=excls->index[a_mol]; j<excls->index[a_mol+1]; j++)
1295                     {
1296                         aj_mol = excls->a[j];
1297                         /* This computation of jla is only correct intra-cg */
1298                         jla = la + aj_mol - a_mol;
1299                         if (jla >= la0 && jla < la1)
1300                         {
1301                             /* This is an intra-cg exclusion. We can skip
1302                              *  the global indexing and distance checking.
1303                              */
1304                             /* Intra-cg exclusions are only required
1305                              * for the home zone.
1306                              */
1307                             if (ic == 0)
1308                             {
1309                                 lexcls->a[n++] = jla;
1310                                 /* Check to avoid double counts */
1311                                 if (jla > la)
1312                                 {
1313                                     count++;
1314                                 }
1315                             }
1316                         }
1317                         else
1318                         {
1319                             /* This is a inter-cg exclusion */
1320                             /* Since exclusions are pair interactions,
1321                              * just like non-bonded interactions,
1322                              * they can be assigned properly up
1323                              * to the DD cutoff (not cutoff_min as
1324                              * for the other bonded interactions).
1325                              */
1326                             if (ga2la_get(ga2la,a_gl+aj_mol-a_mol,&jla,&cell))
1327                             {
1328                                 if (ic == 0 && cell == 0)
1329                                 {
1330                                     lexcls->a[n++] = jla;
1331                                     /* Check to avoid double counts */
1332                                     if (jla > la)
1333                                     {
1334                                         count++;
1335                                     }
1336                                 }
1337                                 else if (jla >= jla0 && jla < jla1 &&
1338                                          (!bRCheck ||
1339                                           dd_dist2(pbc_null,cg_cm,la2lc,la,jla) < rc2))
1340                                 {
1341                                     /* jla > la, since jla0 > la */
1342                                     lexcls->a[n++] = jla;
1343                                     count++;
1344                                 }
1345                             }
1346                         }
1347                     }
1348                 }
1349             }
1350             else
1351             {
1352                 /* There are no inter-cg excls and this cg is self-excluded.
1353                  * These exclusions are only required for zone 0,
1354                  * since other zones do not see themselves.
1355                  */
1356                 if (ic == 0)
1357                 {
1358                     for(la=la0; la<la1; la++)
1359                     {
1360                         lexcls->index[la] = n;
1361                         for(j=la0; j<la1; j++)
1362                         {
1363                             lexcls->a[n++] = j;
1364                         }
1365                     }
1366                     count += ((la1 - la0)*(la1 - la0 - 1))/2;
1367                 }
1368                 else
1369                 {
1370                     /* We don't need exclusions for this cg */
1371                     for(la=la0; la<la1; la++)
1372                     {
1373                         lexcls->index[la] = n;
1374                     }
1375                 }
1376             }
1377         }
1378     }
1379     if (dd->n_intercg_excl == 0)
1380     {
1381         /* There are no exclusions involving non-home charge groups,
1382          * but we need to set the indices for neighborsearching.
1383          */
1384         la0 = dd->cgindex[zones->izone[0].cg1];
1385         for(la=la0; la<lexcls->nr; la++)
1386         {
1387             lexcls->index[la] = n;
1388         }
1389     }
1390     lexcls->index[lexcls->nr] = n;
1391     lexcls->nra = n;
1392     if (dd->n_intercg_excl == 0)
1393     {
1394         /* nr is only used to loop over the exclusions for Ewald and RF,
1395          * so we can set it to the number of home atoms for efficiency.
1396          */
1397         lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1398     }
1399     if (debug)
1400     {
1401         fprintf(debug,"We have %d exclusions, check count %d\n",
1402                 lexcls->nra,count);
1403     }
1404
1405     return count;
1406 }
1407
1408 void dd_make_local_cgs(gmx_domdec_t *dd,t_block *lcgs)
1409 {
1410   lcgs->nr    = dd->ncg_tot;
1411   lcgs->index = dd->cgindex;
1412 }
1413
1414 void dd_make_local_top(FILE *fplog,
1415                        gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1416                        int npbcdim,matrix box,
1417                        rvec cellsize_min,ivec npulse,
1418                        t_forcerec *fr,gmx_vsite_t *vsite,
1419                        gmx_mtop_t *mtop,gmx_localtop_t *ltop)
1420 {
1421     bool bUniqueExcl,bRCheckMB,bRCheck2B,bRCheckExcl;
1422     real rc=-1;
1423     ivec rcheck;
1424     int  d,nexcl;
1425     t_pbc pbc,*pbc_null=NULL;
1426
1427     if (debug)
1428     {
1429         fprintf(debug,"Making local topology\n");
1430     }
1431
1432     dd_make_local_cgs(dd,&ltop->cgs);
1433
1434     bRCheckMB   = FALSE;
1435     bRCheck2B   = FALSE;
1436     bRCheckExcl = FALSE;
1437
1438     if (!dd->reverse_top->bMultiCGmols)
1439     {
1440         /* We don't need checks, assign all interactions with local atoms */
1441
1442         dd->nbonded_local = make_local_bondeds_intracg(dd,mtop->molblock,
1443                                                        &ltop->idef,vsite);
1444     }
1445     else
1446     {
1447         /* We need to check to which cell bondeds should be assigned */
1448         rc = dd_cutoff_twobody(dd);
1449         if (debug)
1450         {
1451             fprintf(debug,"Two-body bonded cut-off distance is %g\n",rc);
1452         }
1453
1454         /* Should we check cg_cm distances when assigning bonded interactions? */
1455         for(d=0; d<DIM; d++)
1456         {
1457             rcheck[d] = FALSE;
1458             /* Only need to check for dimensions where the part of the box
1459              * that is not communicated is smaller than the cut-off.
1460              */
1461             if (d < npbcdim && dd->nc[d] > 1 &&
1462                 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
1463             {
1464                 if (dd->nc[d] == 2)
1465                 {
1466                     rcheck[d] = TRUE;
1467                     bRCheckMB = TRUE;
1468                 }
1469                 /* Check for interactions between two atoms,
1470                  * where we can allow interactions up to the cut-off,
1471                  * instead of up to the smallest cell dimension.
1472                  */
1473                 bRCheck2B = TRUE;
1474             }
1475             if (debug)
1476             {
1477                 fprintf(debug,
1478                         "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1479                         d,cellsize_min[d],d,rcheck[d],bRCheck2B);
1480             }
1481         }
1482         if (dd->reverse_top->bExclRequired)
1483         {
1484             bRCheckExcl = bRCheck2B;
1485         }
1486         else
1487         {
1488             /* If we don't have forces on exclusions,
1489              * we don't care about exclusions being assigned mulitple times.
1490              */
1491             bRCheckExcl = FALSE;
1492         }
1493         if (bRCheckMB || bRCheck2B)
1494         {
1495             make_la2lc(dd);
1496             if (fr->bMolPBC)
1497             {
1498                 set_pbc_dd(&pbc,fr->ePBC,dd,TRUE,box);
1499                 pbc_null = &pbc;
1500             }
1501             else
1502             {
1503                 pbc_null = NULL;
1504             }
1505         }
1506
1507         dd->nbonded_local = make_local_bondeds(dd,zones,mtop->molblock,
1508                                                bRCheckMB,rcheck,bRCheck2B,rc,
1509                                                dd->la2lc,
1510                                                pbc_null,fr->cg_cm,
1511                                                &ltop->idef,vsite);
1512     }
1513
1514     if (dd->reverse_top->ilsort == ilsortNO_FE)
1515     {
1516         ltop->idef.ilsort = ilsortNO_FE;
1517     }
1518     else
1519     {
1520         gmx_sort_ilist_fe(&ltop->idef);
1521     }
1522
1523     nexcl = make_local_exclusions(dd,zones,mtop,bRCheckExcl,
1524                                   rc,dd->la2lc,pbc_null,fr->cg_cm,
1525                                   fr,&ltop->excls);
1526
1527     if (dd->reverse_top->bExclRequired)
1528     {
1529         dd->nbonded_local += nexcl;
1530     }
1531
1532     ltop->atomtypes  = mtop->atomtypes;
1533
1534     /* For an error message only */
1535     dd->reverse_top->err_top_global = mtop;
1536     dd->reverse_top->err_top_local  = ltop;
1537 }
1538
1539 gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
1540 {
1541     gmx_localtop_t *top;
1542     int i;
1543
1544     snew(top,1);
1545
1546     top->idef.ntypes   = top_global->ffparams.ntypes;
1547     top->idef.atnr     = top_global->ffparams.atnr;
1548     top->idef.functype = top_global->ffparams.functype;
1549     top->idef.iparams  = top_global->ffparams.iparams;
1550     top->idef.fudgeQQ  = top_global->ffparams.fudgeQQ;
1551
1552     for(i=0; i<F_NRE; i++)
1553     {
1554         top->idef.il[i].iatoms = NULL;
1555         top->idef.il[i].nalloc = 0;
1556     }
1557     top->idef.ilsort   = ilsortUNKNOWN;
1558
1559     return top;
1560 }
1561
1562 void dd_init_local_state(gmx_domdec_t *dd,
1563                          t_state *state_global,t_state *state_local)
1564 {
1565     int buf[3];
1566
1567     if (DDMASTER(dd))
1568     {
1569         buf[0] = state_global->flags;
1570         buf[1] = state_global->ngtc;
1571         buf[2] = state_global->nnhchains;
1572     }
1573     dd_bcast(dd,3*sizeof(int),buf);
1574
1575     init_state(state_local,0,buf[1],buf[2]);
1576     state_local->flags = buf[0];
1577
1578     /* With Langevin Dynamics we need to make proper storage space
1579      * in the global and local state for the random numbers.
1580      */
1581     if (state_local->flags & (1<<estLD_RNG))
1582     {
1583         if (DDMASTER(dd) && state_global->nrngi > 1)
1584         {
1585             state_global->nrng = dd->nnodes*gmx_rng_n();
1586             srenew(state_global->ld_rng,state_global->nrng);
1587         }
1588         state_local->nrng = gmx_rng_n();
1589         snew(state_local->ld_rng,state_local->nrng);
1590     }
1591     if (state_local->flags & (1<<estLD_RNGI))
1592     {
1593         if (DDMASTER(dd) && state_global->nrngi > 1)
1594         {
1595             state_global->nrngi = dd->nnodes;
1596             srenew(state_global->ld_rngi,state_global->nrngi);
1597         }
1598         snew(state_local->ld_rngi,1);
1599     }
1600 }
1601
1602 static void check_link(t_blocka *link,int cg_gl,int cg_gl_j)
1603 {
1604     int  k,aj;
1605     bool bFound;
1606
1607     bFound = FALSE;
1608     for(k=link->index[cg_gl]; k<link->index[cg_gl+1]; k++)
1609     {
1610         if (link->a[k] == cg_gl_j)
1611         {
1612             bFound = TRUE;
1613         }
1614     }
1615     if (!bFound)
1616     {
1617         /* Add this charge group link */
1618         if (link->index[cg_gl+1]+1 > link->nalloc_a)
1619         {
1620             link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
1621             srenew(link->a,link->nalloc_a);
1622         }
1623         link->a[link->index[cg_gl+1]] = cg_gl_j;
1624         link->index[cg_gl+1]++;
1625     }
1626 }
1627
1628 static int *make_at2cg(t_block *cgs)
1629 {
1630     int *at2cg,cg,a;
1631
1632     snew(at2cg,cgs->index[cgs->nr]);
1633     for(cg=0; cg<cgs->nr; cg++)
1634     {
1635         for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1636         {
1637             at2cg[a] = cg;
1638         }
1639     }
1640
1641     return at2cg;
1642 }
1643
1644 t_blocka *make_charge_group_links(gmx_mtop_t *mtop,gmx_domdec_t *dd,
1645                                   cginfo_mb_t *cginfo_mb)
1646 {
1647     gmx_reverse_top_t *rt;
1648     int  mb,cg_offset,cg,cg_gl,a,aj,i,j,ftype,nral,nlink_mol,mol,ncgi;
1649     gmx_molblock_t *molb;
1650     gmx_moltype_t *molt;
1651     t_block *cgs;
1652     t_blocka *excls;
1653     int *a2c;
1654     gmx_reverse_ilist_t ril;
1655     t_blocka *link;
1656     cginfo_mb_t *cgi_mb;
1657
1658     /* For each charge group make a list of other charge groups
1659      * in the system that a linked to it via bonded interactions
1660      * which are also stored in reverse_top.
1661      */
1662
1663     rt = dd->reverse_top;
1664
1665     snew(link,1);
1666     snew(link->index,ncg_mtop(mtop)+1);
1667     link->nalloc_a = 0;
1668     link->a = NULL;
1669
1670     link->index[0] = 0;
1671     cg_offset = 0;
1672     ncgi = 0;
1673     for(mb=0; mb<mtop->nmolblock; mb++)
1674     {
1675         molb = &mtop->molblock[mb];
1676         if (molb->nmol == 0)
1677         {
1678             continue;
1679         }
1680         molt = &mtop->moltype[molb->type];
1681         cgs   = &molt->cgs;
1682         excls = &molt->excls;
1683         a2c = make_at2cg(cgs);
1684         /* Make a reverse ilist in which the interactions are linked
1685          * to all atoms, not only the first atom as in gmx_reverse_top.
1686          * The constraints are discarded here.
1687          */
1688         make_reverse_ilist(molt,NULL,FALSE,FALSE,TRUE,&ril);
1689
1690         cgi_mb = &cginfo_mb[mb];
1691
1692         for(cg=0; cg<cgs->nr; cg++)
1693         {
1694             cg_gl = cg_offset + cg;
1695             link->index[cg_gl+1] = link->index[cg_gl];
1696             for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1697             {
1698                 i = ril.index[a];
1699                 while (i < ril.index[a+1])
1700                 {
1701                     ftype = ril.il[i++];
1702                     nral = NRAL(ftype);
1703                     /* Skip the ifunc index */
1704                     i++;
1705                     for(j=0; j<nral; j++)
1706                     {
1707                         aj = ril.il[i+j];
1708                         if (a2c[aj] != cg)
1709                         {
1710                             check_link(link,cg_gl,cg_offset+a2c[aj]);
1711                         }
1712                     }
1713                     i += nral_rt(ftype);
1714                 }
1715                 if (rt->bExclRequired)
1716                 {
1717                     /* Exclusions always go both ways */
1718                     for(j=excls->index[a]; j<excls->index[a+1]; j++)
1719                     {
1720                         aj = excls->a[j];
1721                         if (a2c[aj] != cg)
1722                         {
1723                             check_link(link,cg_gl,cg_offset+a2c[aj]);
1724                         }
1725                     }
1726                 }
1727             }
1728             if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
1729             {
1730                 SET_CGINFO_BOND_INTER(cgi_mb[mb].cginfo[cg]);
1731                 ncgi++;
1732             }
1733         }
1734         nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
1735
1736         cg_offset += cgs->nr;
1737
1738         destroy_reverse_ilist(&ril);
1739         sfree(a2c);
1740
1741         if (debug)
1742         {
1743             fprintf(debug,"molecule type '%s' %d cgs has %d cg links through bonded interac.\n",*molt->name,cgs->nr,nlink_mol);
1744         }
1745
1746         if (molb->nmol > 1)
1747         {
1748             /* Copy the data for the rest of the molecules in this block */
1749             link->nalloc_a += (molb->nmol - 1)*nlink_mol;
1750             srenew(link->a,link->nalloc_a);
1751             for(mol=1; mol<molb->nmol; mol++)
1752             {
1753                 for(cg=0; cg<cgs->nr; cg++)
1754                 {
1755                     cg_gl = cg_offset + cg;
1756                     link->index[cg_gl+1] =
1757                         link->index[cg_gl+1-cgs->nr] + nlink_mol;
1758                     for(j=link->index[cg_gl]; j<link->index[cg_gl+1]; j++)
1759                     {
1760                         link->a[j] = link->a[j-nlink_mol] + cgs->nr;
1761                     }
1762                     if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
1763                         cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
1764                     {
1765                         SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
1766                         ncgi++;
1767                     }
1768                 }
1769                 cg_offset += cgs->nr;
1770             }
1771         }
1772     }
1773
1774     if (debug)
1775     {
1776         fprintf(debug,"Of the %d charge groups %d are linked via bonded interactions\n",ncg_mtop(mtop),ncgi);
1777     }
1778
1779     return link;
1780 }
1781
1782 static void bonded_cg_distance_mol(gmx_moltype_t *molt,int *at2cg,
1783                                    bool bBCheck,bool bExcl,rvec *cg_cm,
1784                                    real *r_2b,int *ft2b,int *a2_1,int *a2_2,
1785                                    real *r_mb,int *ftmb,int *am_1,int *am_2)
1786 {
1787     int ftype,nral,i,j,ai,aj,cgi,cgj;
1788     t_ilist *il;
1789     t_blocka *excls;
1790     real r2_2b,r2_mb,rij2;
1791
1792     r2_2b = 0;
1793     r2_mb = 0;
1794     for(ftype=0; ftype<F_NRE; ftype++)
1795     {
1796         if (dd_check_ftype(ftype,bBCheck,FALSE))
1797         {
1798             il = &molt->ilist[ftype];
1799             nral = NRAL(ftype);
1800             if (nral > 1)
1801             {
1802                 for(i=0; i<il->nr; i+=1+nral)
1803                 {
1804                     for(ai=0; ai<nral; ai++)
1805                     {
1806                         cgi = at2cg[il->iatoms[i+1+ai]];
1807                                                 for(aj=0; aj<nral; aj++) {
1808                                                         cgj = at2cg[il->iatoms[i+1+aj]];
1809                                                         if (cgi != cgj)
1810                             {
1811                                                                 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1812                                                                 if (nral == 2 && rij2 > r2_2b)
1813                                 {
1814                                                                         r2_2b = rij2;
1815                                     *ft2b = ftype;
1816                                     *a2_1 = il->iatoms[i+1+ai];
1817                                     *a2_2 = il->iatoms[i+1+aj];
1818                                                                 }
1819                                                                 if (nral >  2 && rij2 > r2_mb)
1820                                 {
1821                                                                         r2_mb = rij2;
1822                                     *ftmb = ftype;
1823                                     *am_1 = il->iatoms[i+1+ai];
1824                                     *am_2 = il->iatoms[i+1+aj];
1825                                                                 }
1826                                                         }
1827                                                 }
1828                                         }
1829                                 }
1830                         }
1831                 }
1832         }
1833         if (bExcl)
1834     {
1835                 excls = &molt->excls;
1836                 for(ai=0; ai<excls->nr; ai++)
1837         {
1838                         cgi = at2cg[ai];
1839                         for(j=excls->index[ai]; j<excls->index[ai+1]; j++) {
1840                                 cgj = at2cg[excls->a[j]];
1841                                 if (cgi != cgj)
1842                 {
1843                                         rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1844                                         if (rij2 > r2_2b)
1845                     {
1846                                                 r2_2b = rij2;
1847                                         }
1848                                 }
1849                         }
1850                 }
1851         }
1852
1853         *r_2b = sqrt(r2_2b);
1854         *r_mb = sqrt(r2_mb);
1855 }
1856
1857 static void get_cgcm_mol(gmx_moltype_t *molt,gmx_ffparams_t *ffparams,
1858                          int ePBC,t_graph *graph,matrix box,
1859                          gmx_vsite_t *vsite,
1860                          rvec *x,rvec *xs,rvec *cg_cm)
1861 {
1862     int n,i;
1863
1864     if (ePBC != epbcNONE)
1865     {
1866         mk_mshift(NULL,graph,ePBC,box,x);
1867
1868         shift_x(graph,box,x,xs);
1869         /* By doing an extra mk_mshift the molecules that are broken
1870          * because they were e.g. imported from another software
1871          * will be made whole again. Such are the healing powers
1872          * of GROMACS.
1873          */
1874         mk_mshift(NULL,graph,ePBC,box,xs);
1875     }
1876     else
1877     {
1878         /* We copy the coordinates so the original coordinates remain
1879          * unchanged, just to be 100% sure that we do not affect
1880          * binary reproducability of simulations.
1881          */
1882         n = molt->cgs.index[molt->cgs.nr];
1883         for(i=0; i<n; i++)
1884         {
1885             copy_rvec(x[i],xs[i]);
1886         }
1887     }
1888
1889     if (vsite)
1890     {
1891         construct_vsites(NULL,vsite,xs,NULL,0.0,NULL,
1892                          ffparams->iparams,molt->ilist,
1893                          epbcNONE,TRUE,NULL,NULL,NULL);
1894     }
1895
1896     calc_cgcm(NULL,0,molt->cgs.nr,&molt->cgs,xs,cg_cm);
1897 }
1898
1899 static int have_vsite_molt(gmx_moltype_t *molt)
1900 {
1901     int  i;
1902     bool bVSite;
1903
1904     bVSite = FALSE;
1905     for(i=0; i<F_NRE; i++)
1906     {
1907         if ((interaction_function[i].flags & IF_VSITE) &&
1908             molt->ilist[i].nr > 0) {
1909             bVSite = TRUE;
1910         }
1911     }
1912
1913     return bVSite;
1914 }
1915
1916 void dd_bonded_cg_distance(FILE *fplog,
1917                            gmx_domdec_t *dd,gmx_mtop_t *mtop,
1918                            t_inputrec *ir,rvec *x,matrix box,
1919                            bool bBCheck,
1920                            real *r_2b,real *r_mb)
1921 {
1922     bool bExclRequired;
1923     int  mb,cg_offset,at_offset,*at2cg,mol;
1924     t_graph graph;
1925     gmx_vsite_t *vsite;
1926     gmx_molblock_t *molb;
1927     gmx_moltype_t *molt;
1928     rvec *xs,*cg_cm;
1929     real rmol_2b,rmol_mb;
1930     int ft2b=-1,a_2b_1=-1,a_2b_2=-1,ftmb=-1,a_mb_1=-1,a_mb_2=-1;
1931     int ftm2b=-1,amol_2b_1=-1,amol_2b_2=-1,ftmmb=-1,amol_mb_1=-1,amol_mb_2=-1;
1932
1933     bExclRequired = IR_EXCL_FORCES(*ir);
1934
1935     /* For gmx_vsite_t everything 0 should work (without pbc) */
1936     snew(vsite,1);
1937
1938     *r_2b = 0;
1939     *r_mb = 0;
1940     cg_offset = 0;
1941     at_offset = 0;
1942     for(mb=0; mb<mtop->nmolblock; mb++)
1943     {
1944         molb = &mtop->molblock[mb];
1945         molt = &mtop->moltype[molb->type];
1946         if (molt->cgs.nr == 1 || molb->nmol == 0)
1947         {
1948             cg_offset += molb->nmol*molt->cgs.nr;
1949             at_offset += molb->nmol*molt->atoms.nr;
1950         }
1951         else
1952         {
1953             if (ir->ePBC != epbcNONE)
1954             {
1955                 mk_graph_ilist(NULL,molt->ilist,0,molt->atoms.nr,FALSE,FALSE,
1956                                &graph);
1957             }
1958
1959             at2cg = make_at2cg(&molt->cgs);
1960             snew(xs,molt->atoms.nr);
1961             snew(cg_cm,molt->cgs.nr);
1962             for(mol=0; mol<molb->nmol; mol++)
1963             {
1964                 get_cgcm_mol(molt,&mtop->ffparams,ir->ePBC,&graph,box,
1965                              have_vsite_molt(molt) ? vsite : NULL,
1966                              x+at_offset,xs,cg_cm);
1967
1968                 bonded_cg_distance_mol(molt,at2cg,bBCheck,bExclRequired,cg_cm,
1969                                        &rmol_2b,&ftm2b,&amol_2b_1,&amol_2b_2,
1970                                        &rmol_mb,&ftmmb,&amol_mb_1,&amol_mb_2);
1971                 if (rmol_2b > *r_2b)
1972                 {
1973                     *r_2b  = rmol_2b;
1974                     ft2b   = ftm2b;
1975                     a_2b_1 = at_offset + amol_2b_1;
1976                     a_2b_2 = at_offset + amol_2b_2;
1977                 }
1978                 if (rmol_mb > *r_mb)
1979                 {
1980                     *r_mb  = rmol_mb;
1981                     ftmb   = ftmmb;
1982                     a_mb_1 = at_offset + amol_mb_1;
1983                     a_mb_2 = at_offset + amol_mb_2;
1984                 }
1985
1986                 cg_offset += molt->cgs.nr;
1987                 at_offset += molt->atoms.nr;
1988             }
1989             sfree(cg_cm);
1990             sfree(xs);
1991             sfree(at2cg);
1992             if (ir->ePBC != epbcNONE)
1993             {
1994                 done_graph(&graph);
1995             }
1996         }
1997     }
1998
1999     sfree(vsite);
2000
2001     if (fplog && (ft2b >= 0 || ftmb >= 0))
2002     {
2003         fprintf(fplog,
2004                 "Initial maximum inter charge-group distances:\n");
2005         if (ft2b >= 0)
2006         {
2007             fprintf(fplog,
2008                     "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2009                     *r_2b,interaction_function[ft2b].longname,
2010                     a_2b_1+1,a_2b_2+1);
2011         }
2012         if (ftmb >= 0)
2013         {
2014             fprintf(fplog,
2015                     "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2016                     *r_mb,interaction_function[ftmb].longname,
2017                     a_mb_1+1,a_mb_2+1);
2018         }
2019     }
2020 }