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