src/mdlib/domdec_top.c

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