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_fbposres(int mol,int a_mol,gmx_molblock_t *molb,
 769                        t_iatom *iatoms,t_idef *idef)
 770 {
 771     int n,a_molb;
 772     t_iparams *ip;
 773
 774     /* This flat-bottom position restraint has not been added yet,
 775      * so it's index is the current number of position restraints.
 776      */
 777     n = idef->il[F_FBPOSRES].nr/2;
 778     if (n+1 > idef->iparams_fbposres_nalloc)
 779     {
 780         idef->iparams_fbposres_nalloc = over_alloc_dd(n+1);
 781         srenew(idef->iparams_fbposres,idef->iparams_fbposres_nalloc);
 782     }
 783     ip = &idef->iparams_fbposres[n];
 784     /* Copy the force constants */
 785     *ip = idef->iparams[iatoms[0]];
 786
 787     /* Get the position restriant coordinats from the molblock */
 788     a_molb = mol*molb->natoms_mol + a_mol;
 789     if (a_molb >= molb->nposres_xA)
 790     {
 791         gmx_incons("Not enough position restraint coordinates");
 792     }
 793     /* Take reference positions from A position of normal posres */
 794     ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 795     ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 796     ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 797
 798     /* Note: no B-type for flat-bottom posres */
 799
 800     /* Set the parameter index for idef->iparams_posre */
 801     iatoms[0] = n;
 802 }
 803
 804 static void add_vsite(gmx_ga2la_t ga2la,int *index,int *rtil,
 805                       int ftype,int nral,
 806                       gmx_bool bHomeA,int a,int a_gl,int a_mol,
 807                       t_iatom *iatoms,
 808                       t_idef *idef,int **vsite_pbc,int *vsite_pbc_nalloc)
 809 {
 810     int  k,ak_gl,vsi,pbc_a_mol;
 811     t_iatom tiatoms[1+MAXATOMLIST],*iatoms_r;
 812     int  j,ftype_r,nral_r;
 813
 814     /* Copy the type */
 815     tiatoms[0] = iatoms[0];
 816
 817     if (bHomeA)
 818     {
 819         /* We know the local index of the first atom */
 820         tiatoms[1] = a;
 821     }
 822     else
 823     {
 824         /* Convert later in make_local_vsites */
 825         tiatoms[1] = -a_gl - 1;
 826     }
 827
 828     for(k=2; k<1+nral; k++)
 829     {
 830         ak_gl = a_gl + iatoms[k] - a_mol;
 831         if (!ga2la_get_home(ga2la,ak_gl,&tiatoms[k]))
 832         {
 833             /* Copy the global index, convert later in make_local_vsites */
 834             tiatoms[k] = -(ak_gl + 1);
 835         }
 836     }
 837
 838     /* Add this interaction to the local topology */
 839     add_ifunc(nral,tiatoms,&idef->il[ftype]);
 840     if (vsite_pbc)
 841     {
 842         vsi = idef->il[ftype].nr/(1+nral) - 1;
 843         if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
 844         {
 845             vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
 846             srenew(vsite_pbc[ftype-F_VSITE2],vsite_pbc_nalloc[ftype-F_VSITE2]);
 847         }
 848         if (bHomeA)
 849         {
 850             pbc_a_mol = iatoms[1+nral+1];
 851             if (pbc_a_mol < 0)
 852             {
 853                 /* The pbc flag is one of the following two options:
 854                  * -2: vsite and all constructing atoms are within the same cg, no pbc
 855                  * -1: vsite and its first constructing atom are in the same cg, do pbc
 856                  */
 857                 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
 858             }
 859             else
 860             {
 861                 /* Set the pbc atom for this vsite so we can make its pbc
 862                  * identical to the rest of the atoms in its charge group.
 863                  * Since the order of the atoms does not change within a charge
 864                  * group, we do not need the global to local atom index.
 865                  */
 866                 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
 867             }
 868         }
 869         else
 870         {
 871             /* This vsite is non-home (required for recursion),
 872              * and therefore there is no charge group to match pbc with.
 873              * But we always turn on full_pbc to assure that higher order
 874              * recursion works correctly.
 875              */
 876             vsite_pbc[ftype-F_VSITE2][vsi] = -1;
 877         }
 878     }
 879
 880     if (iatoms[1+nral])
 881     {
 882         /* Check for recursion */
 883         for(k=2; k<1+nral; k++)
 884         {
 885             if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
 886             {
 887                 /* This construction atoms is a vsite and not a home atom */
 888                 if (gmx_debug_at)
 889                 {
 890                     fprintf(debug,"Constructing atom %d of vsite atom %d is a vsite and non-home\n",iatoms[k]+1,a_mol+1);
 891                 }
 892                 /* Find the vsite construction */
 893
 894                 /* Check all interactions assigned to this atom */
 895                 j = index[iatoms[k]];
 896                 while (j < index[iatoms[k]+1])
 897                 {
 898                     ftype_r = rtil[j++];
 899                     nral_r = NRAL(ftype_r);
 900                     if (interaction_function[ftype_r].flags & IF_VSITE)
 901                     {
 902                         /* Add this vsite (recursion) */
 903                         add_vsite(ga2la,index,rtil,ftype_r,nral_r,
 904                                   FALSE,-1,a_gl+iatoms[k]-iatoms[1],iatoms[k],
 905                                   rtil+j,idef,vsite_pbc,vsite_pbc_nalloc);
 906                         j += 1 + nral_r + 2;
 907                     }
 908                     else
 909                     {
 910                         j += 1 + nral_r;
 911                     }
 912                 }
 913             }
 914         }
 915     }
 916 }
 917
 918 static void make_la2lc(gmx_domdec_t *dd)
 919 {
 920     int *cgindex,*la2lc,cg,a;
 921
 922     cgindex = dd->cgindex;
 923
 924     if (dd->nat_tot > dd->la2lc_nalloc)
 925     {
 926         dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
 927         snew(dd->la2lc,dd->la2lc_nalloc);
 928     }
 929     la2lc = dd->la2lc;
 930
 931     /* Make the local atom to local cg index */
 932     for(cg=0; cg<dd->ncg_tot; cg++)
 933     {
 934         for(a=cgindex[cg]; a<cgindex[cg+1]; a++)
 935         {
 936             la2lc[a] = cg;
 937         }
 938     }
 939 }
 940
 941 static real dd_dist2(t_pbc *pbc_null,rvec *cg_cm,const int *la2lc,int i,int j)
 942 {
 943     rvec dx;
 944
 945     if (pbc_null)
 946     {
 947         pbc_dx_aiuc(pbc_null,cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
 948     }
 949     else
 950     {
 951         rvec_sub(cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
 952     }
 953
 954     return norm2(dx);
 955 }
 956
 957 static int make_local_bondeds(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
 958                               gmx_molblock_t *molb,
 959                               gmx_bool bRCheckMB,ivec rcheck,gmx_bool bRCheck2B,
 960                               real rc,
 961                               int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
 962                               t_idef *idef,gmx_vsite_t *vsite)
 963 {
 964     int nzone,nizone,ic,la0,la1,i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,d,k;
 965     int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
 966     t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
 967     gmx_bool bBCheck,bUse,bLocal;
 968     real rc2;
 969     ivec k_zero,k_plus;
 970     gmx_ga2la_t ga2la;
 971     int  a_loc;
 972     int  kc;
 973     gmx_domdec_ns_ranges_t *izone;
 974     gmx_reverse_top_t *rt;
 975     gmx_molblock_ind_t *mbi;
 976     int nbonded_local;
 977
 978     nzone  = zones->n;
 979     nizone = zones->nizone;
 980     izone  = zones->izone;
 981
 982     rc2 = rc*rc;
 983
 984     if (vsite && vsite->n_intercg_vsite > 0)
 985     {
 986         vsite_pbc        = vsite->vsite_pbc_loc;
 987         vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
 988     }
 989     else
 990     {
 991         vsite_pbc        = NULL;
 992         vsite_pbc_nalloc = NULL;
 993     }
 994
 995     rt = dd->reverse_top;
 996
 997     bBCheck = rt->bBCheck;
 998
 999     /* Clear the counts */
1000     for(ftype=0; ftype<F_NRE; ftype++)
1001     {
1002         idef->il[ftype].nr = 0;
1003     }
1004     nbonded_local = 0;
1005
1006     mbi = rt->mbi;
1007
1008     ga2la = dd->ga2la;
1009
1010     for(ic=0; ic<nzone; ic++)
1011     {
1012         la0 = dd->cgindex[zones->cg_range[ic]];
1013         la1 = dd->cgindex[zones->cg_range[ic+1]];
1014         for(i=la0; i<la1; i++)
1015         {
1016             /* Get the global atom number */
1017             i_gl = dd->gatindex[i];
1018             global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
1019             /* Check all interactions assigned to this atom */
1020             index = rt->ril_mt[mt].index;
1021             rtil  = rt->ril_mt[mt].il;
1022             j = index[i_mol];
1023             while (j < index[i_mol+1])
1024             {
1025                 ftype  = rtil[j++];
1026                 iatoms = rtil + j;
1027                 nral = NRAL(ftype);
1028                 if (interaction_function[ftype].flags & IF_VSITE)
1029                 {
1030                     /* The vsite construction goes where the vsite itself is */
1031                     if (ic == 0)
1032                     {
1033                         add_vsite(dd->ga2la,index,rtil,ftype,nral,
1034                                   TRUE,i,i_gl,i_mol,
1035                                   iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1036                     }
1037                     j += 1 + nral + 2;
1038                 }
1039                 else
1040                 {
1041                     /* Copy the type */
1042                     tiatoms[0] = iatoms[0];
1043
1044                     if (nral == 1)
1045                     {
1046                         /* Assign single-body interactions to the home zone */
1047                         if (ic == 0)
1048                         {
1049                             bUse = TRUE;
1050                             tiatoms[1] = i;
1051                             if (ftype == F_POSRES)
1052                             {
1053                                 add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1054                             }
1055                             else if (ftype == F_FBPOSRES)
1056                             {
1057                                 add_fbposres(mol,i_mol,&molb[mb],tiatoms,idef);
1058                             }
1059                         }
1060                         else
1061                         {
1062                             bUse = FALSE;
1063                         }
1064                     }
1065                     else if (nral == 2)
1066                     {
1067                         /* This is a two-body interaction, we can assign
1068                          * analogous to the non-bonded assignments.
1069                          */
1070                         if (!ga2la_get(ga2la,i_gl+iatoms[2]-i_mol,&a_loc,&kc))
1071                         {
1072                             bUse = FALSE;
1073                         }
1074                         else
1075                         {
1076                             if (kc >= nzone)
1077                             {
1078                                 kc -= nzone;
1079                             }
1080                             /* Check zone interaction assignments */
1081                             bUse = ((ic < nizone && ic <= kc &&
1082                                      izone[ic].j0 <= kc && kc < izone[ic].j1) ||
1083                                     (kc < nizone && ic >  kc &&
1084                                      izone[kc].j0 <= ic && ic < izone[kc].j1));
1085                             if (bUse)
1086                             {
1087                                 tiatoms[1] = i;
1088                                 tiatoms[2] = a_loc;
1089                                 /* If necessary check the cgcm distance */
1090                                 if (bRCheck2B &&
1091                                     dd_dist2(pbc_null,cg_cm,la2lc,
1092                                              tiatoms[1],tiatoms[2]) >= rc2)
1093                                 {
1094                                     bUse = FALSE;
1095                                 }
1096                             }
1097                         }
1098                     }
1099                     else
1100                     {
1101                         /* Assign this multi-body bonded interaction to
1102                          * the local node if we have all the atoms involved
1103                          * (local or communicated) and the minimum zone shift
1104                          * in each dimension is zero, for dimensions
1105                          * with 2 DD cells an extra check may be necessary.
1106                          */
1107                         bUse = TRUE;
1108                         clear_ivec(k_zero);
1109                         clear_ivec(k_plus);
1110                         for(k=1; k<=nral && bUse; k++)
1111                         {
1112                             bLocal = ga2la_get(ga2la,i_gl+iatoms[k]-i_mol,
1113                                                &a_loc,&kc);
1114                             if (!bLocal || kc >= zones->n)
1115                             {
1116                                 /* We do not have this atom of this interaction
1117                                  * locally, or it comes from more than one cell
1118                                  * away.
1119                                  */
1120                                 bUse = FALSE;
1121                             }
1122                             else
1123                             {
1124                                 tiatoms[k] = a_loc;
1125                                 for(d=0; d<DIM; d++)
1126                                 {
1127                                     if (zones->shift[kc][d] == 0)
1128                                     {
1129                                         k_zero[d] = k;
1130                                     }
1131                                     else
1132                                     {
1133                                         k_plus[d] = k;
1134                                     }
1135                                 }
1136                             }
1137                         }
1138                         bUse = (bUse &&
1139                                 k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1140                         if (bRCheckMB)
1141                         {
1142                             for(d=0; (d<DIM && bUse); d++)
1143                             {
1144                                 /* Check if the cg_cm distance falls within
1145                                  * the cut-off to avoid possible multiple
1146                                  * assignments of bonded interactions.
1147                                  */
1148                                 if (rcheck[d] &&
1149                                     k_plus[d] &&
1150                                     dd_dist2(pbc_null,cg_cm,la2lc,
1151                                              tiatoms[k_zero[d]],tiatoms[k_plus[d]]) >= rc2)
1152                                 {
1153                                     bUse = FALSE;
1154                                 }
1155                             }
1156                         }
1157                     }
1158                     if (bUse)
1159                     {
1160                         /* Add this interaction to the local topology */
1161                         add_ifunc(nral,tiatoms,&idef->il[ftype]);
1162                         /* Sum so we can check in global_stat
1163                          * if we have everything.
1164                          */
1165                         if (bBCheck ||
1166                             !(interaction_function[ftype].flags & IF_LIMZERO))
1167                         {
1168                             nbonded_local++;
1169                         }
1170                     }
1171                     j += 1 + nral;
1172                 }
1173             }
1174         }
1175     }
1176
1177     return nbonded_local;
1178 }
1179
1180 static int make_local_bondeds_intracg(gmx_domdec_t *dd,gmx_molblock_t *molb,
1181                                       t_idef *idef,gmx_vsite_t *vsite)
1182 {
1183     int i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,k;
1184     int *index,*rtil,**vsite_pbc,*vsite_pbc_nalloc;
1185     t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
1186     gmx_reverse_top_t *rt;
1187     gmx_molblock_ind_t *mbi;
1188     int nbonded_local;
1189
1190     if (vsite && vsite->n_intercg_vsite > 0)
1191     {
1192         vsite_pbc        = vsite->vsite_pbc_loc;
1193         vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1194     }
1195     else
1196     {
1197         vsite_pbc        = NULL;
1198         vsite_pbc_nalloc = NULL;
1199     }
1200
1201     /* Clear the counts */
1202     for(ftype=0; ftype<F_NRE; ftype++)
1203     {
1204         idef->il[ftype].nr = 0;
1205     }
1206     nbonded_local = 0;
1207
1208     rt = dd->reverse_top;
1209
1210     if (rt->ril_mt_tot_size == 0)
1211     {
1212         /* There are no interactions to assign */
1213         return nbonded_local;
1214     }
1215
1216     mbi = rt->mbi;
1217
1218     for(i=0; i<dd->nat_home; i++)
1219     {
1220         /* Get the global atom number */
1221         i_gl = dd->gatindex[i];
1222         global_atomnr_to_moltype_ind(mbi,i_gl,&mb,&mt,&mol,&i_mol);
1223         /* Check all interactions assigned to this atom */
1224         index = rt->ril_mt[mt].index;
1225         rtil  = rt->ril_mt[mt].il;
1226         /* Check all interactions assigned to this atom */
1227         j = index[i_mol];
1228         while (j < index[i_mol+1])
1229         {
1230             ftype  = rtil[j++];
1231             iatoms = rtil + j;
1232             nral = NRAL(ftype);
1233             if (interaction_function[ftype].flags & IF_VSITE)
1234             {
1235                 /* The vsite construction goes where the vsite itself is */
1236                 add_vsite(dd->ga2la,index,rtil,ftype,nral,
1237                           TRUE,i,i_gl,i_mol,
1238                           iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1239                 j += 1 + nral + 2;
1240             }
1241             else
1242             {
1243                 /* Copy the type */
1244                 tiatoms[0] = iatoms[0];
1245                 tiatoms[1] = i;
1246                 for(k=2; k<=nral; k++)
1247                 {
1248                     tiatoms[k] = i + iatoms[k] - iatoms[1];
1249                 }
1250                 if (ftype == F_POSRES)
1251                 {
1252                     add_posres(mol,i_mol,&molb[mb],tiatoms,idef);
1253                 }
1254                 else if (ftype == F_FBPOSRES)
1255                 {
1256                     add_fbposres(mol,i_mol,&molb[mb],tiatoms,idef);
1257                 }
1258                 /* Add this interaction to the local topology */
1259                 add_ifunc(nral,tiatoms,&idef->il[ftype]);
1260                 /* Sum so we can check in global_stat if we have everything */
1261                 nbonded_local++;
1262                 j += 1 + nral;
1263             }
1264         }
1265     }
1266
1267     return nbonded_local;
1268 }
1269
1270 static int make_local_exclusions(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1271                                  gmx_mtop_t *mtop,
1272                                  gmx_bool bRCheck,real rc,
1273                                  int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1274                                  t_forcerec *fr,
1275                                  t_blocka *lexcls)
1276 {
1277     int  nizone,n,count,ic,jla0,jla1,jla;
1278     int  cg,la0,la1,la,a_gl,mb,mt,mol,a_mol,j,aj_mol;
1279     t_blocka *excls;
1280     gmx_ga2la_t ga2la;
1281     int  a_loc;
1282     int  cell;
1283     gmx_molblock_ind_t *mbi;
1284     real rc2;
1285
1286     /* Since for RF and PME we need to loop over the exclusions
1287      * we should store each exclusion only once. This is done
1288      * using the same zone scheme as used for neighbor searching.
1289      * The exclusions involving non-home atoms are stored only
1290      * one way: atom j is in the excl list of i only for j > i,
1291      * where i and j are local atom numbers.
1292      */
1293
1294     lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1295     if (lexcls->nr+1 > lexcls->nalloc_index)
1296     {
1297         lexcls->nalloc_index = over_alloc_dd(lexcls->nr)+1;
1298         srenew(lexcls->index,lexcls->nalloc_index);
1299     }
1300
1301     mbi = dd->reverse_top->mbi;
1302
1303     ga2la = dd->ga2la;
1304
1305     rc2 = rc*rc;
1306
1307     if (dd->n_intercg_excl)
1308     {
1309         nizone = zones->nizone;
1310     }
1311     else
1312     {
1313         nizone = 1;
1314     }
1315     n = 0;
1316     count = 0;
1317     for(ic=0; ic<nizone; ic++)
1318     {
1319         jla0 = dd->cgindex[zones->izone[ic].jcg0];
1320         jla1 = dd->cgindex[zones->izone[ic].jcg1];
1321         for(cg=zones->cg_range[ic]; cg<zones->cg_range[ic+1]; cg++)
1322         {
1323             /* Here we assume the number of exclusions in one charge group
1324              * is never larger than 1000.
1325              */
1326             if (n+1000 > lexcls->nalloc_a)
1327             {
1328                 lexcls->nalloc_a = over_alloc_large(n+1000);
1329                 srenew(lexcls->a,lexcls->nalloc_a);
1330             }
1331             la0 = dd->cgindex[cg];
1332             la1 = dd->cgindex[cg+1];
1333             if (GET_CGINFO_EXCL_INTER(fr->cginfo[cg]) ||
1334                 !GET_CGINFO_EXCL_INTRA(fr->cginfo[cg]))
1335             {
1336                 /* Copy the exclusions from the global top */
1337                 for(la=la0; la<la1; la++) {
1338                     lexcls->index[la] = n;
1339                     a_gl = dd->gatindex[la];
1340                     global_atomnr_to_moltype_ind(mbi,a_gl,&mb,&mt,&mol,&a_mol);
1341                     excls = &mtop->moltype[mt].excls;
1342                     for(j=excls->index[a_mol]; j<excls->index[a_mol+1]; j++)
1343                     {
1344                         aj_mol = excls->a[j];
1345                         /* This computation of jla is only correct intra-cg */
1346                         jla = la + aj_mol - a_mol;
1347                         if (jla >= la0 && jla < la1)
1348                         {
1349                             /* This is an intra-cg exclusion. We can skip
1350                              *  the global indexing and distance checking.
1351                              */
1352                             /* Intra-cg exclusions are only required
1353                              * for the home zone.
1354                              */
1355                             if (ic == 0)
1356                             {
1357                                 lexcls->a[n++] = jla;
1358                                 /* Check to avoid double counts */
1359                                 if (jla > la)
1360                                 {
1361                                     count++;
1362                                 }
1363                             }
1364                         }
1365                         else
1366                         {
1367                             /* This is a inter-cg exclusion */
1368                             /* Since exclusions are pair interactions,
1369                              * just like non-bonded interactions,
1370                              * they can be assigned properly up
1371                              * to the DD cutoff (not cutoff_min as
1372                              * for the other bonded interactions).
1373                              */
1374                             if (ga2la_get(ga2la,a_gl+aj_mol-a_mol,&jla,&cell))
1375                             {
1376                                 if (ic == 0 && cell == 0)
1377                                 {
1378                                     lexcls->a[n++] = jla;
1379                                     /* Check to avoid double counts */
1380                                     if (jla > la)
1381                                     {
1382                                         count++;
1383                                     }
1384                                 }
1385                                 else if (jla >= jla0 && jla < jla1 &&
1386                                          (!bRCheck ||
1387                                           dd_dist2(pbc_null,cg_cm,la2lc,la,jla) < rc2))
1388                                 {
1389                                     /* jla > la, since jla0 > la */
1390                                     lexcls->a[n++] = jla;
1391                                     count++;
1392                                 }
1393                             }
1394                         }
1395                     }
1396                 }
1397             }
1398             else
1399             {
1400                 /* There are no inter-cg excls and this cg is self-excluded.
1401                  * These exclusions are only required for zone 0,
1402                  * since other zones do not see themselves.
1403                  */
1404                 if (ic == 0)
1405                 {
1406                     for(la=la0; la<la1; la++)
1407                     {
1408                         lexcls->index[la] = n;
1409                         for(j=la0; j<la1; j++)
1410                         {
1411                             lexcls->a[n++] = j;
1412                         }
1413                     }
1414                     count += ((la1 - la0)*(la1 - la0 - 1))/2;
1415                 }
1416                 else
1417                 {
1418                     /* We don't need exclusions for this cg */
1419                     for(la=la0; la<la1; la++)
1420                     {
1421                         lexcls->index[la] = n;
1422                     }
1423                 }
1424             }
1425         }
1426     }
1427     if (dd->n_intercg_excl == 0)
1428     {
1429         /* There are no exclusions involving non-home charge groups,
1430          * but we need to set the indices for neighborsearching.
1431          */
1432         la0 = dd->cgindex[zones->izone[0].cg1];
1433         for(la=la0; la<lexcls->nr; la++)
1434         {
1435             lexcls->index[la] = n;
1436         }
1437     }
1438     lexcls->index[lexcls->nr] = n;
1439     lexcls->nra = n;
1440     if (dd->n_intercg_excl == 0)
1441     {
1442         /* nr is only used to loop over the exclusions for Ewald and RF,
1443          * so we can set it to the number of home atoms for efficiency.
1444          */
1445         lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1446     }
1447     if (debug)
1448     {
1449         fprintf(debug,"We have %d exclusions, check count %d\n",
1450                 lexcls->nra,count);
1451     }
1452
1453     return count;
1454 }
1455
1456 void dd_make_local_cgs(gmx_domdec_t *dd,t_block *lcgs)
1457 {
1458   lcgs->nr    = dd->ncg_tot;
1459   lcgs->index = dd->cgindex;
1460 }
1461
1462 void dd_make_local_top(FILE *fplog,
1463                        gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1464                        int npbcdim,matrix box,
1465                        rvec cellsize_min,ivec npulse,
1466                        t_forcerec *fr,gmx_vsite_t *vsite,
1467                        gmx_mtop_t *mtop,gmx_localtop_t *ltop)
1468 {
1469     gmx_bool bUniqueExcl,bRCheckMB,bRCheck2B,bRCheckExcl;
1470     real rc=-1;
1471     ivec rcheck;
1472     int  d,nexcl;
1473     t_pbc pbc,*pbc_null=NULL;
1474
1475     if (debug)
1476     {
1477         fprintf(debug,"Making local topology\n");
1478     }
1479
1480     dd_make_local_cgs(dd,&ltop->cgs);
1481
1482     bRCheckMB   = FALSE;
1483     bRCheck2B   = FALSE;
1484     bRCheckExcl = FALSE;
1485
1486     if (!dd->reverse_top->bMultiCGmols)
1487     {
1488         /* We don't need checks, assign all interactions with local atoms */
1489
1490         dd->nbonded_local = make_local_bondeds_intracg(dd,mtop->molblock,
1491                                                        &ltop->idef,vsite);
1492     }
1493     else
1494     {
1495         /* We need to check to which cell bondeds should be assigned */
1496         rc = dd_cutoff_twobody(dd);
1497         if (debug)
1498         {
1499             fprintf(debug,"Two-body bonded cut-off distance is %g\n",rc);
1500         }
1501
1502         /* Should we check cg_cm distances when assigning bonded interactions? */
1503         for(d=0; d<DIM; d++)
1504         {
1505             rcheck[d] = FALSE;
1506             /* Only need to check for dimensions where the part of the box
1507              * that is not communicated is smaller than the cut-off.
1508              */
1509             if (d < npbcdim && dd->nc[d] > 1 &&
1510                 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
1511             {
1512                 if (dd->nc[d] == 2)
1513                 {
1514                     rcheck[d] = TRUE;
1515                     bRCheckMB = TRUE;
1516                 }
1517                 /* Check for interactions between two atoms,
1518                  * where we can allow interactions up to the cut-off,
1519                  * instead of up to the smallest cell dimension.
1520                  */
1521                 bRCheck2B = TRUE;
1522             }
1523             if (debug)
1524             {
1525                 fprintf(debug,
1526                         "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1527                         d,cellsize_min[d],d,rcheck[d],bRCheck2B);
1528             }
1529         }
1530         if (dd->reverse_top->bExclRequired)
1531         {
1532             bRCheckExcl = bRCheck2B;
1533         }
1534         else
1535         {
1536             /* If we don't have forces on exclusions,
1537              * we don't care about exclusions being assigned mulitple times.
1538              */
1539             bRCheckExcl = FALSE;
1540         }
1541         if (bRCheckMB || bRCheck2B)
1542         {
1543             make_la2lc(dd);
1544             if (fr->bMolPBC)
1545             {
1546                 set_pbc_dd(&pbc,fr->ePBC,dd,TRUE,box);
1547                 pbc_null = &pbc;
1548             }
1549             else
1550             {
1551                 pbc_null = NULL;
1552             }
1553         }
1554
1555         dd->nbonded_local = make_local_bondeds(dd,zones,mtop->molblock,
1556                                                bRCheckMB,rcheck,bRCheck2B,rc,
1557                                                dd->la2lc,
1558                                                pbc_null,fr->cg_cm,
1559                                                &ltop->idef,vsite);
1560     }
1561
1562     /* The ilist is not sorted yet,
1563      * we can only do this when we have the charge arrays.
1564      */
1565     ltop->idef.ilsort = ilsortUNKNOWN;
1566
1567     nexcl = make_local_exclusions(dd,zones,mtop,bRCheckExcl,
1568                                   rc,dd->la2lc,pbc_null,fr->cg_cm,
1569                                   fr,&ltop->excls);
1570
1571     if (dd->reverse_top->bExclRequired)
1572     {
1573         dd->nbonded_local += nexcl;
1574     }
1575
1576     ltop->atomtypes  = mtop->atomtypes;
1577
1578     /* For an error message only */
1579     dd->reverse_top->err_top_global = mtop;
1580     dd->reverse_top->err_top_local  = ltop;
1581 }
1582
1583 void dd_sort_local_top(gmx_domdec_t *dd,t_mdatoms *mdatoms,
1584                        gmx_localtop_t *ltop)
1585 {
1586     if (dd->reverse_top->ilsort == ilsortNO_FE)
1587     {
1588         ltop->idef.ilsort = ilsortNO_FE;
1589     }
1590     else
1591     {
1592         gmx_sort_ilist_fe(&ltop->idef,mdatoms->chargeA,mdatoms->chargeB);
1593     }
1594 }
1595
1596 gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
1597 {
1598     gmx_localtop_t *top;
1599     int i;
1600
1601     snew(top,1);
1602
1603     top->idef.ntypes   = top_global->ffparams.ntypes;
1604     top->idef.atnr     = top_global->ffparams.atnr;
1605     top->idef.functype = top_global->ffparams.functype;
1606     top->idef.iparams  = top_global->ffparams.iparams;
1607     top->idef.fudgeQQ  = top_global->ffparams.fudgeQQ;
1608     top->idef.cmap_grid= top_global->ffparams.cmap_grid;
1609
1610     for(i=0; i<F_NRE; i++)
1611     {
1612         top->idef.il[i].iatoms = NULL;
1613         top->idef.il[i].nalloc = 0;
1614     }
1615     top->idef.ilsort   = ilsortUNKNOWN;
1616
1617     return top;
1618 }
1619
1620 void dd_init_local_state(gmx_domdec_t *dd,
1621                          t_state *state_global,t_state *state_local)
1622 {
1623     int buf[NITEM_DD_INIT_LOCAL_STATE];
1624
1625     if (DDMASTER(dd))
1626     {
1627         buf[0] = state_global->flags;
1628         buf[1] = state_global->ngtc;
1629         buf[2] = state_global->nnhpres;
1630         buf[3] = state_global->nhchainlength;
1631         buf[4] = state_global->dfhist.nlambda;
1632     }
1633     dd_bcast(dd,NITEM_DD_INIT_LOCAL_STATE*sizeof(int),buf);
1634
1635     init_state(state_local,0,buf[1],buf[2],buf[3],buf[4]);
1636     state_local->flags = buf[0];
1637
1638     /* With Langevin Dynamics we need to make proper storage space
1639      * in the global and local state for the random numbers.
1640      */
1641     if (state_local->flags & (1<<estLD_RNG))
1642     {
1643         if (DDMASTER(dd) && state_global->nrngi > 1)
1644         {
1645             state_global->nrng = dd->nnodes*gmx_rng_n();
1646             srenew(state_global->ld_rng,state_global->nrng);
1647         }
1648         state_local->nrng = gmx_rng_n();
1649         snew(state_local->ld_rng,state_local->nrng);
1650     }
1651     if (state_local->flags & (1<<estLD_RNGI))
1652     {
1653         if (DDMASTER(dd) && state_global->nrngi > 1)
1654         {
1655             state_global->nrngi = dd->nnodes;
1656             srenew(state_global->ld_rngi,state_global->nrngi);
1657         }
1658         snew(state_local->ld_rngi,1);
1659     }
1660 }
1661
1662 static void check_link(t_blocka *link,int cg_gl,int cg_gl_j)
1663 {
1664     int  k,aj;
1665     gmx_bool bFound;
1666
1667     bFound = FALSE;
1668     for(k=link->index[cg_gl]; k<link->index[cg_gl+1]; k++)
1669     {
1670         if (link->a[k] == cg_gl_j)
1671         {
1672             bFound = TRUE;
1673         }
1674     }
1675     if (!bFound)
1676     {
1677         /* Add this charge group link */
1678         if (link->index[cg_gl+1]+1 > link->nalloc_a)
1679         {
1680             link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
1681             srenew(link->a,link->nalloc_a);
1682         }
1683         link->a[link->index[cg_gl+1]] = cg_gl_j;
1684         link->index[cg_gl+1]++;
1685     }
1686 }
1687
1688 static int *make_at2cg(t_block *cgs)
1689 {
1690     int *at2cg,cg,a;
1691
1692     snew(at2cg,cgs->index[cgs->nr]);
1693     for(cg=0; cg<cgs->nr; cg++)
1694     {
1695         for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1696         {
1697             at2cg[a] = cg;
1698         }
1699     }
1700
1701     return at2cg;
1702 }
1703
1704 t_blocka *make_charge_group_links(gmx_mtop_t *mtop,gmx_domdec_t *dd,
1705                                   cginfo_mb_t *cginfo_mb)
1706 {
1707     gmx_reverse_top_t *rt;
1708     int  mb,cg_offset,cg,cg_gl,a,aj,i,j,ftype,nral,nlink_mol,mol,ncgi;
1709     gmx_molblock_t *molb;
1710     gmx_moltype_t *molt;
1711     t_block *cgs;
1712     t_blocka *excls;
1713     int *a2c;
1714     gmx_reverse_ilist_t ril;
1715     t_blocka *link;
1716     cginfo_mb_t *cgi_mb;
1717
1718     /* For each charge group make a list of other charge groups
1719      * in the system that a linked to it via bonded interactions
1720      * which are also stored in reverse_top.
1721      */
1722
1723     rt = dd->reverse_top;
1724
1725     snew(link,1);
1726     snew(link->index,ncg_mtop(mtop)+1);
1727     link->nalloc_a = 0;
1728     link->a = NULL;
1729
1730     link->index[0] = 0;
1731     cg_offset = 0;
1732     ncgi = 0;
1733     for(mb=0; mb<mtop->nmolblock; mb++)
1734     {
1735         molb = &mtop->molblock[mb];
1736         if (molb->nmol == 0)
1737         {
1738             continue;
1739         }
1740         molt = &mtop->moltype[molb->type];
1741         cgs   = &molt->cgs;
1742         excls = &molt->excls;
1743         a2c = make_at2cg(cgs);
1744         /* Make a reverse ilist in which the interactions are linked
1745          * to all atoms, not only the first atom as in gmx_reverse_top.
1746          * The constraints are discarded here.
1747          */
1748         make_reverse_ilist(molt,NULL,FALSE,FALSE,TRUE,&ril);
1749
1750         cgi_mb = &cginfo_mb[mb];
1751
1752         for(cg=0; cg<cgs->nr; cg++)
1753         {
1754             cg_gl = cg_offset + cg;
1755             link->index[cg_gl+1] = link->index[cg_gl];
1756             for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
1757             {
1758                 i = ril.index[a];
1759                 while (i < ril.index[a+1])
1760                 {
1761                     ftype = ril.il[i++];
1762                     nral = NRAL(ftype);
1763                     /* Skip the ifunc index */
1764                     i++;
1765                     for(j=0; j<nral; j++)
1766                     {
1767                         aj = ril.il[i+j];
1768                         if (a2c[aj] != cg)
1769                         {
1770                             check_link(link,cg_gl,cg_offset+a2c[aj]);
1771                         }
1772                     }
1773                     i += nral_rt(ftype);
1774                 }
1775                 if (rt->bExclRequired)
1776                 {
1777                     /* Exclusions always go both ways */
1778                     for(j=excls->index[a]; j<excls->index[a+1]; j++)
1779                     {
1780                         aj = excls->a[j];
1781                         if (a2c[aj] != cg)
1782                         {
1783                             check_link(link,cg_gl,cg_offset+a2c[aj]);
1784                         }
1785                     }
1786                 }
1787             }
1788             if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
1789             {
1790                 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
1791                 ncgi++;
1792             }
1793         }
1794         nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
1795
1796         cg_offset += cgs->nr;
1797
1798         destroy_reverse_ilist(&ril);
1799         sfree(a2c);
1800
1801         if (debug)
1802         {
1803             fprintf(debug,"molecule type '%s' %d cgs has %d cg links through bonded interac.\n",*molt->name,cgs->nr,nlink_mol);
1804         }
1805
1806         if (molb->nmol > 1)
1807         {
1808             /* Copy the data for the rest of the molecules in this block */
1809             link->nalloc_a += (molb->nmol - 1)*nlink_mol;
1810             srenew(link->a,link->nalloc_a);
1811             for(mol=1; mol<molb->nmol; mol++)
1812             {
1813                 for(cg=0; cg<cgs->nr; cg++)
1814                 {
1815                     cg_gl = cg_offset + cg;
1816                     link->index[cg_gl+1] =
1817                         link->index[cg_gl+1-cgs->nr] + nlink_mol;
1818                     for(j=link->index[cg_gl]; j<link->index[cg_gl+1]; j++)
1819                     {
1820                         link->a[j] = link->a[j-nlink_mol] + cgs->nr;
1821                     }
1822                     if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
1823                         cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
1824                     {
1825                         SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
1826                         ncgi++;
1827                     }
1828                 }
1829                 cg_offset += cgs->nr;
1830             }
1831         }
1832     }
1833
1834     if (debug)
1835     {
1836         fprintf(debug,"Of the %d charge groups %d are linked via bonded interactions\n",ncg_mtop(mtop),ncgi);
1837     }
1838
1839     return link;
1840 }
1841
1842 static void bonded_cg_distance_mol(gmx_moltype_t *molt,int *at2cg,
1843                                    gmx_bool bBCheck,gmx_bool bExcl,rvec *cg_cm,
1844                                    real *r_2b,int *ft2b,int *a2_1,int *a2_2,
1845                                    real *r_mb,int *ftmb,int *am_1,int *am_2)
1846 {
1847     int ftype,nral,i,j,ai,aj,cgi,cgj;
1848     t_ilist *il;
1849     t_blocka *excls;
1850     real r2_2b,r2_mb,rij2;
1851
1852     r2_2b = 0;
1853     r2_mb = 0;
1854     for(ftype=0; ftype<F_NRE; ftype++)
1855     {
1856         if (dd_check_ftype(ftype,bBCheck,FALSE))
1857         {
1858             il = &molt->ilist[ftype];
1859             nral = NRAL(ftype);
1860             if (nral > 1)
1861             {
1862                 for(i=0; i<il->nr; i+=1+nral)
1863                 {
1864                     for(ai=0; ai<nral; ai++)
1865                     {
1866                         cgi = at2cg[il->iatoms[i+1+ai]];
1867                                                 for(aj=0; aj<nral; aj++) {
1868                                                         cgj = at2cg[il->iatoms[i+1+aj]];
1869                                                         if (cgi != cgj)
1870                             {
1871                                                                 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1872                                                                 if (nral == 2 && rij2 > r2_2b)
1873                                 {
1874                                                                         r2_2b = rij2;
1875                                     *ft2b = ftype;
1876                                     *a2_1 = il->iatoms[i+1+ai];
1877                                     *a2_2 = il->iatoms[i+1+aj];
1878                                                                 }
1879                                                                 if (nral >  2 && rij2 > r2_mb)
1880                                 {
1881                                                                         r2_mb = rij2;
1882                                     *ftmb = ftype;
1883                                     *am_1 = il->iatoms[i+1+ai];
1884                                     *am_2 = il->iatoms[i+1+aj];
1885                                                                 }
1886                                                         }
1887                                                 }
1888                                         }
1889                                 }
1890                         }
1891                 }
1892         }
1893         if (bExcl)
1894     {
1895                 excls = &molt->excls;
1896                 for(ai=0; ai<excls->nr; ai++)
1897         {
1898                         cgi = at2cg[ai];
1899                         for(j=excls->index[ai]; j<excls->index[ai+1]; j++) {
1900                                 cgj = at2cg[excls->a[j]];
1901                                 if (cgi != cgj)
1902                 {
1903                                         rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
1904                                         if (rij2 > r2_2b)
1905                     {
1906                                                 r2_2b = rij2;
1907                                         }
1908                                 }
1909                         }
1910                 }
1911         }
1912
1913         *r_2b = sqrt(r2_2b);
1914         *r_mb = sqrt(r2_mb);
1915 }
1916
1917 static void get_cgcm_mol(gmx_moltype_t *molt,gmx_ffparams_t *ffparams,
1918                          int ePBC,t_graph *graph,matrix box,
1919                          gmx_vsite_t *vsite,
1920                          rvec *x,rvec *xs,rvec *cg_cm)
1921 {
1922     int n,i;
1923
1924     if (ePBC != epbcNONE)
1925     {
1926         mk_mshift(NULL,graph,ePBC,box,x);
1927
1928         shift_x(graph,box,x,xs);
1929         /* By doing an extra mk_mshift the molecules that are broken
1930          * because they were e.g. imported from another software
1931          * will be made whole again. Such are the healing powers
1932          * of GROMACS.
1933          */
1934         mk_mshift(NULL,graph,ePBC,box,xs);
1935     }
1936     else
1937     {
1938         /* We copy the coordinates so the original coordinates remain
1939          * unchanged, just to be 100% sure that we do not affect
1940          * binary reproducibility of simulations.
1941          */
1942         n = molt->cgs.index[molt->cgs.nr];
1943         for(i=0; i<n; i++)
1944         {
1945             copy_rvec(x[i],xs[i]);
1946         }
1947     }
1948
1949     if (vsite)
1950     {
1951         construct_vsites(NULL,vsite,xs,NULL,0.0,NULL,
1952                          ffparams->iparams,molt->ilist,
1953                          epbcNONE,TRUE,NULL,NULL,NULL);
1954     }
1955
1956     calc_cgcm(NULL,0,molt->cgs.nr,&molt->cgs,xs,cg_cm);
1957 }
1958
1959 static int have_vsite_molt(gmx_moltype_t *molt)
1960 {
1961     int  i;
1962     gmx_bool bVSite;
1963
1964     bVSite = FALSE;
1965     for(i=0; i<F_NRE; i++)
1966     {
1967         if ((interaction_function[i].flags & IF_VSITE) &&
1968             molt->ilist[i].nr > 0) {
1969             bVSite = TRUE;
1970         }
1971     }
1972
1973     return bVSite;
1974 }
1975
1976 void dd_bonded_cg_distance(FILE *fplog,
1977                            gmx_domdec_t *dd,gmx_mtop_t *mtop,
1978                            t_inputrec *ir,rvec *x,matrix box,
1979                            gmx_bool bBCheck,
1980                            real *r_2b,real *r_mb)
1981 {
1982     gmx_bool bExclRequired;
1983     int  mb,cg_offset,at_offset,*at2cg,mol;
1984     t_graph graph;
1985     gmx_vsite_t *vsite;
1986     gmx_molblock_t *molb;
1987     gmx_moltype_t *molt;
1988     rvec *xs,*cg_cm;
1989     real rmol_2b,rmol_mb;
1990     int ft2b=-1,a_2b_1=-1,a_2b_2=-1,ftmb=-1,a_mb_1=-1,a_mb_2=-1;
1991     int ftm2b=-1,amol_2b_1=-1,amol_2b_2=-1,ftmmb=-1,amol_mb_1=-1,amol_mb_2=-1;
1992
1993     bExclRequired = IR_EXCL_FORCES(*ir);
1994
1995     /* For gmx_vsite_t everything 0 should work (without pbc) */
1996     snew(vsite,1);
1997
1998     *r_2b = 0;
1999     *r_mb = 0;
2000     cg_offset = 0;
2001     at_offset = 0;
2002     for(mb=0; mb<mtop->nmolblock; mb++)
2003     {
2004         molb = &mtop->molblock[mb];
2005         molt = &mtop->moltype[molb->type];
2006         if (molt->cgs.nr == 1 || molb->nmol == 0)
2007         {
2008             cg_offset += molb->nmol*molt->cgs.nr;
2009             at_offset += molb->nmol*molt->atoms.nr;
2010         }
2011         else
2012         {
2013             if (ir->ePBC != epbcNONE)
2014             {
2015                 mk_graph_ilist(NULL,molt->ilist,0,molt->atoms.nr,FALSE,FALSE,
2016                                &graph);
2017             }
2018
2019             at2cg = make_at2cg(&molt->cgs);
2020             snew(xs,molt->atoms.nr);
2021             snew(cg_cm,molt->cgs.nr);
2022             for(mol=0; mol<molb->nmol; mol++)
2023             {
2024                 get_cgcm_mol(molt,&mtop->ffparams,ir->ePBC,&graph,box,
2025                              have_vsite_molt(molt) ? vsite : NULL,
2026                              x+at_offset,xs,cg_cm);
2027
2028                 bonded_cg_distance_mol(molt,at2cg,bBCheck,bExclRequired,cg_cm,
2029                                        &rmol_2b,&ftm2b,&amol_2b_1,&amol_2b_2,
2030                                        &rmol_mb,&ftmmb,&amol_mb_1,&amol_mb_2);
2031                 if (rmol_2b > *r_2b)
2032                 {
2033                     *r_2b  = rmol_2b;
2034                     ft2b   = ftm2b;
2035                     a_2b_1 = at_offset + amol_2b_1;
2036                     a_2b_2 = at_offset + amol_2b_2;
2037                 }
2038                 if (rmol_mb > *r_mb)
2039                 {
2040                     *r_mb  = rmol_mb;
2041                     ftmb   = ftmmb;
2042                     a_mb_1 = at_offset + amol_mb_1;
2043                     a_mb_2 = at_offset + amol_mb_2;
2044                 }
2045
2046                 cg_offset += molt->cgs.nr;
2047                 at_offset += molt->atoms.nr;
2048             }
2049             sfree(cg_cm);
2050             sfree(xs);
2051             sfree(at2cg);
2052             if (ir->ePBC != epbcNONE)
2053             {
2054                 done_graph(&graph);
2055             }
2056         }
2057     }
2058
2059     sfree(vsite);
2060
2061     if (fplog && (ft2b >= 0 || ftmb >= 0))
2062     {
2063         fprintf(fplog,
2064                 "Initial maximum inter charge-group distances:\n");
2065         if (ft2b >= 0)
2066         {
2067             fprintf(fplog,
2068                     "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2069                     *r_2b,interaction_function[ft2b].longname,
2070                     a_2b_1+1,a_2b_2+1);
2071         }
2072         if (ftmb >= 0)
2073         {
2074             fprintf(fplog,
2075                     "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2076                     *r_mb,interaction_function[ftmb].longname,
2077                     a_mb_1+1,a_mb_2+1);
2078         }
2079     }
2080 }