src/mdlib/domdec_top.c

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