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