src/gromacs/mdlib/domdec_top.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * This file is part of Gromacs        Copyright (c) 1991-2008
   5  * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
   6  *
   7  * This program is free software; you can redistribute it and/or
   8  * modify it under the terms of the GNU General Public License
   9  * as published by the Free Software Foundation; either version 2
  10  * of the License, or (at your option) any later version.
  11  *
  12  * To help us fund GROMACS development, we humbly ask that you cite
  13  * the research papers on the package. Check out http://www.gromacs.org
  14  *
  15  * And Hey:
  16  * Gnomes, ROck Monsters And Chili Sauce
  17  */
  18
  19 #ifdef HAVE_CONFIG_H
  20 #include <config.h>
  21 #endif
  22
  23 #include <string.h>
  24 #include "typedefs.h"
  25 #include "smalloc.h"
  26 #include "domdec.h"
  27 #include "domdec_network.h"
  28 #include "names.h"
  29 #include "network.h"
  30 #include "vec.h"
  31 #include "pbc.h"
  32 #include "chargegroup.h"
  33 #include "gmx_random.h"
  34 #include "topsort.h"
  35 #include "mtop_util.h"
  36 #include "mshift.h"
  37 #include "vsite.h"
  38 #include "gmx_ga2la.h"
  39 #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_fbposres(int mol,int a_mol,const gmx_molblock_t *molb,
 831                          t_iatom *iatoms,const t_iparams *ip_in,
 832                          t_idef *idef)
 833 {
 834     int n,a_molb;
 835     t_iparams *ip;
 836
 837     /* This flat-bottom position restraint has not been added yet,
 838      * so it's index is the current number of position restraints.
 839      */
 840     n = idef->il[F_FBPOSRES].nr/2;
 841     if (n+1 > idef->iparams_fbposres_nalloc)
 842     {
 843         idef->iparams_fbposres_nalloc = over_alloc_dd(n+1);
 844         srenew(idef->iparams_fbposres,idef->iparams_fbposres_nalloc);
 845     }
 846     ip = &idef->iparams_fbposres[n];
 847     /* Copy the force constants */
 848     *ip = ip_in[iatoms[0]];
 849
 850     /* Get the position restriant coordinats from the molblock */
 851     a_molb = mol*molb->natoms_mol + a_mol;
 852     if (a_molb >= molb->nposres_xA)
 853     {
 854         gmx_incons("Not enough position restraint coordinates");
 855     }
 856     /* Take reference positions from A position of normal posres */
 857     ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
 858     ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
 859     ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
 860
 861     /* Note: no B-type for flat-bottom posres */
 862
 863     /* Set the parameter index for idef->iparams_posre */
 864     iatoms[0] = n;
 865 }
 866
 867 static void add_vsite(gmx_ga2la_t ga2la,int *index,int *rtil,
 868                       int ftype,int nral,
 869                       gmx_bool bHomeA,int a,int a_gl,int a_mol,
 870                       t_iatom *iatoms,
 871                       t_idef *idef,int **vsite_pbc,int *vsite_pbc_nalloc)
 872 {
 873     int  k,ak_gl,vsi,pbc_a_mol;
 874     t_iatom tiatoms[1+MAXATOMLIST],*iatoms_r;
 875     int  j,ftype_r,nral_r;
 876
 877     /* Copy the type */
 878     tiatoms[0] = iatoms[0];
 879
 880     if (bHomeA)
 881     {
 882         /* We know the local index of the first atom */
 883         tiatoms[1] = a;
 884     }
 885     else
 886     {
 887         /* Convert later in make_local_vsites */
 888         tiatoms[1] = -a_gl - 1;
 889     }
 890
 891     for(k=2; k<1+nral; k++)
 892     {
 893         ak_gl = a_gl + iatoms[k] - a_mol;
 894         if (!ga2la_get_home(ga2la,ak_gl,&tiatoms[k]))
 895         {
 896             /* Copy the global index, convert later in make_local_vsites */
 897             tiatoms[k] = -(ak_gl + 1);
 898         }
 899     }
 900
 901     /* Add this interaction to the local topology */
 902     add_ifunc(nral,tiatoms,&idef->il[ftype]);
 903     if (vsite_pbc)
 904     {
 905         vsi = idef->il[ftype].nr/(1+nral) - 1;
 906         if (vsi >= vsite_pbc_nalloc[ftype-F_VSITE2])
 907         {
 908             vsite_pbc_nalloc[ftype-F_VSITE2] = over_alloc_large(vsi+1);
 909             srenew(vsite_pbc[ftype-F_VSITE2],vsite_pbc_nalloc[ftype-F_VSITE2]);
 910         }
 911         if (bHomeA)
 912         {
 913             pbc_a_mol = iatoms[1+nral+1];
 914             if (pbc_a_mol < 0)
 915             {
 916                 /* The pbc flag is one of the following two options:
 917                  * -2: vsite and all constructing atoms are within the same cg, no pbc
 918                  * -1: vsite and its first constructing atom are in the same cg, do pbc
 919                  */
 920                 vsite_pbc[ftype-F_VSITE2][vsi] = pbc_a_mol;
 921             }
 922             else
 923             {
 924                 /* Set the pbc atom for this vsite so we can make its pbc
 925                  * identical to the rest of the atoms in its charge group.
 926                  * Since the order of the atoms does not change within a charge
 927                  * group, we do not need the global to local atom index.
 928                  */
 929                 vsite_pbc[ftype-F_VSITE2][vsi] = a + pbc_a_mol - iatoms[1];
 930             }
 931         }
 932         else
 933         {
 934             /* This vsite is non-home (required for recursion),
 935              * and therefore there is no charge group to match pbc with.
 936              * But we always turn on full_pbc to assure that higher order
 937              * recursion works correctly.
 938              */
 939             vsite_pbc[ftype-F_VSITE2][vsi] = -1;
 940         }
 941     }
 942
 943     if (iatoms[1+nral])
 944     {
 945         /* Check for recursion */
 946         for(k=2; k<1+nral; k++)
 947         {
 948             if ((iatoms[1+nral] & (2<<k)) && (tiatoms[k] < 0))
 949             {
 950                 /* This construction atoms is a vsite and not a home atom */
 951                 if (gmx_debug_at)
 952                 {
 953                     fprintf(debug,"Constructing atom %d of vsite atom %d is a vsite and non-home\n",iatoms[k]+1,a_mol+1);
 954                 }
 955                 /* Find the vsite construction */
 956
 957                 /* Check all interactions assigned to this atom */
 958                 j = index[iatoms[k]];
 959                 while (j < index[iatoms[k]+1])
 960                 {
 961                     ftype_r = rtil[j++];
 962                     nral_r = NRAL(ftype_r);
 963                     if (interaction_function[ftype_r].flags & IF_VSITE)
 964                     {
 965                         /* Add this vsite (recursion) */
 966                         add_vsite(ga2la,index,rtil,ftype_r,nral_r,
 967                                   FALSE,-1,a_gl+iatoms[k]-iatoms[1],iatoms[k],
 968                                   rtil+j,idef,vsite_pbc,vsite_pbc_nalloc);
 969                         j += 1 + nral_r + 2;
 970                     }
 971                     else
 972                     {
 973                         j += 1 + nral_r;
 974                     }
 975                 }
 976             }
 977         }
 978     }
 979 }
 980
 981 static void make_la2lc(gmx_domdec_t *dd)
 982 {
 983     int *cgindex,*la2lc,cg,a;
 984
 985     cgindex = dd->cgindex;
 986
 987     if (dd->nat_tot > dd->la2lc_nalloc)
 988     {
 989         dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
 990         snew(dd->la2lc,dd->la2lc_nalloc);
 991     }
 992     la2lc = dd->la2lc;
 993
 994     /* Make the local atom to local cg index */
 995     for(cg=0; cg<dd->ncg_tot; cg++)
 996     {
 997         for(a=cgindex[cg]; a<cgindex[cg+1]; a++)
 998         {
 999             la2lc[a] = cg;
1000         }
1001     }
1002 }
1003
1004 static real dd_dist2(t_pbc *pbc_null,rvec *cg_cm,const int *la2lc,int i,int j)
1005 {
1006     rvec dx;
1007
1008     if (pbc_null)
1009     {
1010         pbc_dx_aiuc(pbc_null,cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
1011     }
1012     else
1013     {
1014         rvec_sub(cg_cm[la2lc[i]],cg_cm[la2lc[j]],dx);
1015     }
1016
1017     return norm2(dx);
1018 }
1019
1020 /* Append the nsrc t_blocka block structures in src to *dest */
1021 static void combine_blocka(t_blocka *dest,const t_blocka *src,int nsrc)
1022 {
1023     int ni,na,s,i;
1024
1025     ni = src[nsrc-1].nr;
1026     na = 0;
1027     for(s=0; s<nsrc; s++)
1028     {
1029         na += src[s].nra;
1030     }
1031     if (ni + 1 > dest->nalloc_index)
1032     {
1033         dest->nalloc_index = over_alloc_large(ni+1);
1034         srenew(dest->index,dest->nalloc_index);
1035     }
1036     if (dest->nra + na > dest->nalloc_a)
1037     {
1038         dest->nalloc_a = over_alloc_large(dest->nra+na);
1039         srenew(dest->a,dest->nalloc_a);
1040     }
1041     for(s=0; s<nsrc; s++)
1042     {
1043         for(i=dest->nr+1; i<src[s].nr+1; i++)
1044         {
1045             dest->index[i] = dest->nra + src[s].index[i];
1046         }
1047         for(i=0; i<src[s].nra; i++)
1048         {
1049             dest->a[dest->nra+i] = src[s].a[i];
1050         }
1051         dest->nr   = src[s].nr;
1052         dest->nra += src[s].nra;
1053     }
1054 }
1055
1056 /* Append the nsrc t_idef structures in src to *dest,
1057  * virtual sites need special attention, as pbc info differs per vsite.
1058  */
1059 static void combine_idef(t_idef *dest,const t_idef *src,int nsrc,
1060                          gmx_vsite_t *vsite,int ***vsite_pbc_t)
1061 {
1062     int ftype,n,s,i;
1063     t_ilist *ild;
1064     const t_ilist *ils;
1065     gmx_bool vpbc;
1066     int nral1=0,ftv=0;
1067
1068     for(ftype=0; ftype<F_NRE; ftype++)
1069     {
1070         n = 0;
1071         for(s=0; s<nsrc; s++)
1072         {
1073             n += src[s].il[ftype].nr;
1074         }
1075         if (n > 0)
1076         {
1077             ild = &dest->il[ftype];
1078
1079             if (ild->nr + n > ild->nalloc)
1080             {
1081                 ild->nalloc = over_alloc_large(ild->nr+n);
1082                 srenew(ild->iatoms,ild->nalloc);
1083             }
1084
1085             vpbc = ((interaction_function[ftype].flags & IF_VSITE) &&
1086                     vsite->vsite_pbc_loc != NULL);
1087             if (vpbc)
1088             {
1089                 nral1 = 1 + NRAL(ftype);
1090                 ftv = ftype - F_VSITE2;
1091                 if ((ild->nr + n)/nral1 > vsite->vsite_pbc_loc_nalloc[ftv])
1092                 {
1093                     vsite->vsite_pbc_loc_nalloc[ftv] =
1094                         over_alloc_large((ild->nr + n)/nral1);
1095                     srenew(vsite->vsite_pbc_loc[ftv],
1096                            vsite->vsite_pbc_loc_nalloc[ftv]);
1097                 }
1098             }
1099
1100             for(s=0; s<nsrc; s++)
1101             {
1102                 ils = &src[s].il[ftype];
1103                 for(i=0; i<ils->nr; i++)
1104                 {
1105                     ild->iatoms[ild->nr+i] = ils->iatoms[i];
1106                 }
1107                 if (vpbc)
1108                 {
1109                     for(i=0; i<ils->nr; i+=nral1)
1110                     {
1111                         vsite->vsite_pbc_loc[ftv][(ild->nr+i)/nral1] =
1112                             vsite_pbc_t[s][ftv][i/nral1];
1113                     }
1114                 }
1115
1116                 ild->nr += ils->nr;
1117             }
1118         }
1119     }
1120
1121     /* Position restraints need an additional treatment */
1122     if (dest->il[F_POSRES].nr > 0)
1123     {
1124         n = dest->il[F_POSRES].nr/2;
1125         if (n > dest->iparams_posres_nalloc)
1126         {
1127             dest->iparams_posres_nalloc = over_alloc_large(n);
1128             srenew(dest->iparams_posres,dest->iparams_posres_nalloc);
1129         }
1130         /* Set n to the number of original position restraints in dest */
1131         for(s=0; s<nsrc; s++)
1132         {
1133             n -= src[s].il[F_POSRES].nr/2;
1134         }
1135         for(s=0; s<nsrc; s++)
1136         {
1137             for(i=0; i<src[s].il[F_POSRES].nr/2; i++)
1138             {
1139                 /* Correct the index into iparams_posres */
1140                 dest->il[F_POSRES].iatoms[n*2] = n;
1141                 /* Copy the position restraint force parameters */
1142                 dest->iparams_posres[n] = src[s].iparams_posres[i];
1143                 n++;
1144             }
1145         }
1146     }
1147 }
1148
1149 /* This function looks up and assigns bonded interactions for zone iz.
1150  * With thread parallelizing each thread acts on a different atom range:
1151  * at_start to at_end.
1152  */
1153 static int make_bondeds_zone(gmx_domdec_t *dd,
1154                              const gmx_domdec_zones_t *zones,
1155                              const gmx_molblock_t *molb,
1156                              gmx_bool bRCheckMB,ivec rcheck,gmx_bool bRCheck2B,
1157                              real rc2,
1158                              int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1159                              const t_iparams *ip_in,
1160                              t_idef *idef,gmx_vsite_t *vsite,
1161                              int **vsite_pbc,
1162                              int *vsite_pbc_nalloc,
1163                              int iz,int nzone,
1164                              int at_start,int at_end)
1165 {
1166     int i,i_gl,mb,mt,mol,i_mol,j,ftype,nral,d,k;
1167     int *index,*rtil;
1168     t_iatom *iatoms,tiatoms[1+MAXATOMLIST];
1169     gmx_bool bBCheck,bUse,bLocal;
1170     ivec k_zero,k_plus;
1171     gmx_ga2la_t ga2la;
1172     int  a_loc;
1173     int  kz;
1174     int  nizone;
1175     const gmx_domdec_ns_ranges_t *izone;
1176     gmx_reverse_top_t *rt;
1177     int nbonded_local;
1178
1179     nizone = zones->nizone;
1180     izone  = zones->izone;
1181
1182     rt = dd->reverse_top;
1183
1184     bBCheck = rt->bBCheck;
1185
1186     nbonded_local = 0;
1187
1188     ga2la = dd->ga2la;
1189
1190     for(i=at_start; i<at_end; i++)
1191     {
1192         /* Get the global atom number */
1193         i_gl = dd->gatindex[i];
1194         global_atomnr_to_moltype_ind(rt,i_gl,&mb,&mt,&mol,&i_mol);
1195         /* Check all interactions assigned to this atom */
1196         index = rt->ril_mt[mt].index;
1197         rtil  = rt->ril_mt[mt].il;
1198         j = index[i_mol];
1199         while (j < index[i_mol+1])
1200         {
1201             ftype  = rtil[j++];
1202             iatoms = rtil + j;
1203             nral = NRAL(ftype);
1204             if (ftype == F_SETTLE)
1205             {
1206                 /* Settles are only in the reverse top when they
1207                  * operate within a charge group. So we can assign
1208                  * them without checks. We do this only for performance
1209                  * reasons; it could be handled by the code below.
1210                  */
1211                 if (iz == 0)
1212                 {
1213                     /* Home zone: add this settle to the local topology */
1214                     tiatoms[0] = iatoms[0];
1215                     tiatoms[1] = i;
1216                     tiatoms[2] = i + iatoms[2] - iatoms[1];
1217                     tiatoms[3] = i + iatoms[3] - iatoms[1];
1218                     add_ifunc(nral,tiatoms,&idef->il[ftype]);
1219                     nbonded_local++;
1220                 }
1221                 j += 1 + nral;
1222             }
1223             else if (interaction_function[ftype].flags & IF_VSITE)
1224             {
1225                 /* The vsite construction goes where the vsite itself is */
1226                 if (iz == 0)
1227                 {
1228                     add_vsite(dd->ga2la,index,rtil,ftype,nral,
1229                               TRUE,i,i_gl,i_mol,
1230                               iatoms,idef,vsite_pbc,vsite_pbc_nalloc);
1231                 }
1232                 j += 1 + nral + 2;
1233             }
1234             else
1235             {
1236                 /* Copy the type */
1237                 tiatoms[0] = iatoms[0];
1238
1239                 if (nral == 1)
1240                 {
1241                     /* Assign single-body interactions to the home zone */
1242                     if (iz == 0)
1243                     {
1244                         bUse = TRUE;
1245                             tiatoms[1] = i;
1246                             if (ftype == F_POSRES)
1247                             {
1248                                 add_posres(mol,i_mol,&molb[mb],tiatoms,ip_in,
1249                                            idef);
1250                             }
1251                             else if (ftype == F_FBPOSRES)
1252                             {
1253                                 add_fbposres(mol,i_mol,&molb[mb],tiatoms,ip_in,
1254                                              idef);
1255                             }
1256                     }
1257                     else
1258                     {
1259                         bUse = FALSE;
1260                     }
1261                 }
1262                 else if (nral == 2)
1263                 {
1264                     /* This is a two-body interaction, we can assign
1265                      * analogous to the non-bonded assignments.
1266                      */
1267                     if (!ga2la_get(ga2la,i_gl+iatoms[2]-i_mol,&a_loc,&kz))
1268                     {
1269                         bUse = FALSE;
1270                     }
1271                     else
1272                     {
1273                         if (kz >= nzone)
1274                         {
1275                             kz -= nzone;
1276                         }
1277                         /* Check zone interaction assignments */
1278                         bUse = ((iz < nizone && iz <= kz &&
1279                                  izone[iz].j0 <= kz && kz < izone[iz].j1) ||
1280                                 (kz < nizone && iz >  kz &&
1281                                  izone[kz].j0 <= iz && iz < izone[kz].j1));
1282                         if (bUse)
1283                         {
1284                             tiatoms[1] = i;
1285                             tiatoms[2] = a_loc;
1286                             /* If necessary check the cgcm distance */
1287                             if (bRCheck2B &&
1288                                 dd_dist2(pbc_null,cg_cm,la2lc,
1289                                          tiatoms[1],tiatoms[2]) >= rc2)
1290                             {
1291                                 bUse = FALSE;
1292                             }
1293                         }
1294                     }
1295                 }
1296                 else
1297                 {
1298                     /* Assign this multi-body bonded interaction to
1299                      * the local node if we have all the atoms involved
1300                      * (local or communicated) and the minimum zone shift
1301                      * in each dimension is zero, for dimensions
1302                      * with 2 DD cells an extra check may be necessary.
1303                      */
1304                     bUse = TRUE;
1305                     clear_ivec(k_zero);
1306                     clear_ivec(k_plus);
1307                     for(k=1; k<=nral && bUse; k++)
1308                     {
1309                         bLocal = ga2la_get(ga2la,i_gl+iatoms[k]-i_mol,
1310                                            &a_loc,&kz);
1311                         if (!bLocal || kz >= zones->n)
1312                         {
1313                             /* We do not have this atom of this interaction
1314                              * locally, or it comes from more than one cell
1315                              * away.
1316                              */
1317                             bUse = FALSE;
1318                         }
1319                         else
1320                         {
1321                             tiatoms[k] = a_loc;
1322                             for(d=0; d<DIM; d++)
1323                             {
1324                                 if (zones->shift[kz][d] == 0)
1325                                 {
1326                                     k_zero[d] = k;
1327                                 }
1328                                 else
1329                                 {
1330                                     k_plus[d] = k;
1331                                 }
1332                             }
1333                         }
1334                     }
1335                     bUse = (bUse &&
1336                             k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
1337                     if (bRCheckMB)
1338                     {
1339                         for(d=0; (d<DIM && bUse); d++)
1340                         {
1341                             /* Check if the cg_cm distance falls within
1342                              * the cut-off to avoid possible multiple
1343                              * assignments of bonded interactions.
1344                              */
1345                             if (rcheck[d] &&
1346                                 k_plus[d] &&
1347                                 dd_dist2(pbc_null,cg_cm,la2lc,
1348                                          tiatoms[k_zero[d]],tiatoms[k_plus[d]]) >= rc2)
1349                             {
1350                                 bUse = FALSE;
1351                             }
1352                         }
1353                     }
1354                 }
1355                 if (bUse)
1356                 {
1357                     /* Add this interaction to the local topology */
1358                     add_ifunc(nral,tiatoms,&idef->il[ftype]);
1359                     /* Sum so we can check in global_stat
1360                      * if we have everything.
1361                      */
1362                     if (bBCheck ||
1363                         !(interaction_function[ftype].flags & IF_LIMZERO))
1364                     {
1365                         nbonded_local++;
1366                     }
1367                 }
1368                 j += 1 + nral;
1369             }
1370         }
1371     }
1372
1373     return nbonded_local;
1374 }
1375
1376 static void set_no_exclusions_zone(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1377                                    int iz,t_blocka *lexcls)
1378 {
1379     int  a0,a1,a;
1380
1381     a0 = dd->cgindex[zones->cg_range[iz]];
1382     a1 = dd->cgindex[zones->cg_range[iz+1]];
1383
1384     for(a=a0+1; a<a1+1; a++)
1385     {
1386         lexcls->index[a] = lexcls->nra;
1387     }
1388 }
1389
1390 static int make_exclusions_zone(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1391                                 const gmx_moltype_t *moltype,
1392                                 gmx_bool bRCheck,real rc2,
1393                                 int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1394                                 const int *cginfo,
1395                                 t_blocka *lexcls,
1396                                 int iz,
1397                                 int cg_start,int cg_end)
1398 {
1399     int  nizone,n,count,jla0,jla1,jla;
1400     int  cg,la0,la1,la,a_gl,mb,mt,mol,a_mol,j,aj_mol;
1401     const t_blocka *excls;
1402     gmx_ga2la_t ga2la;
1403     int  a_loc;
1404     int  cell;
1405
1406     ga2la = dd->ga2la;
1407
1408     jla0 = dd->cgindex[zones->izone[iz].jcg0];
1409     jla1 = dd->cgindex[zones->izone[iz].jcg1];
1410
1411     /* We set the end index, but note that we might not start at zero here */
1412     lexcls->nr = dd->cgindex[cg_end];
1413
1414     n = lexcls->nra;
1415     count = 0;
1416     for(cg=cg_start; cg<cg_end; cg++)
1417     {
1418         /* Here we assume the number of exclusions in one charge group
1419          * is never larger than 1000.
1420          */
1421         if (n+1000 > lexcls->nalloc_a)
1422         {
1423             lexcls->nalloc_a = over_alloc_large(n+1000);
1424             srenew(lexcls->a,lexcls->nalloc_a);
1425         }
1426         la0 = dd->cgindex[cg];
1427         la1 = dd->cgindex[cg+1];
1428         if (GET_CGINFO_EXCL_INTER(cginfo[cg]) ||
1429             !GET_CGINFO_EXCL_INTRA(cginfo[cg]))
1430         {
1431             /* Copy the exclusions from the global top */
1432             for(la=la0; la<la1; la++) {
1433                 lexcls->index[la] = n;
1434                 a_gl = dd->gatindex[la];
1435                 global_atomnr_to_moltype_ind(dd->reverse_top,a_gl,&mb,&mt,&mol,&a_mol);
1436                 excls = &moltype[mt].excls;
1437                 for(j=excls->index[a_mol]; j<excls->index[a_mol+1]; j++)
1438                 {
1439                     aj_mol = excls->a[j];
1440                     /* This computation of jla is only correct intra-cg */
1441                     jla = la + aj_mol - a_mol;
1442                     if (jla >= la0 && jla < la1)
1443                     {
1444                         /* This is an intra-cg exclusion. We can skip
1445                          *  the global indexing and distance checking.
1446                          */
1447                         /* Intra-cg exclusions are only required
1448                          * for the home zone.
1449                          */
1450                         if (iz == 0)
1451                         {
1452                             lexcls->a[n++] = jla;
1453                             /* Check to avoid double counts */
1454                             if (jla > la)
1455                             {
1456                                 count++;
1457                             }
1458                         }
1459                     }
1460                     else
1461                     {
1462                         /* This is a inter-cg exclusion */
1463                         /* Since exclusions are pair interactions,
1464                          * just like non-bonded interactions,
1465                          * they can be assigned properly up
1466                          * to the DD cutoff (not cutoff_min as
1467                          * for the other bonded interactions).
1468                          */
1469                         if (ga2la_get(ga2la,a_gl+aj_mol-a_mol,&jla,&cell))
1470                         {
1471                             if (iz == 0 && cell == 0)
1472                             {
1473                                 lexcls->a[n++] = jla;
1474                                 /* Check to avoid double counts */
1475                                 if (jla > la)
1476                                 {
1477                                     count++;
1478                                 }
1479                             }
1480                             else if (jla >= jla0 && jla < jla1 &&
1481                                      (!bRCheck ||
1482                                       dd_dist2(pbc_null,cg_cm,la2lc,la,jla) < rc2))
1483                             {
1484                                 /* jla > la, since jla0 > la */
1485                                 lexcls->a[n++] = jla;
1486                                 count++;
1487                             }
1488                         }
1489                     }
1490                 }
1491             }
1492         }
1493         else
1494         {
1495             /* There are no inter-cg excls and this cg is self-excluded.
1496              * These exclusions are only required for zone 0,
1497              * since other zones do not see themselves.
1498              */
1499             if (iz == 0)
1500             {
1501                 for(la=la0; la<la1; la++)
1502                 {
1503                     lexcls->index[la] = n;
1504                     for(j=la0; j<la1; j++)
1505                     {
1506                         lexcls->a[n++] = j;
1507                     }
1508                 }
1509                 count += ((la1 - la0)*(la1 - la0 - 1))/2;
1510             }
1511             else
1512             {
1513                 /* We don't need exclusions for this cg */
1514                 for(la=la0; la<la1; la++)
1515                 {
1516                     lexcls->index[la] = n;
1517                 }
1518             }
1519         }
1520     }
1521
1522     lexcls->index[lexcls->nr] = n;
1523     lexcls->nra = n;
1524
1525     return count;
1526 }
1527
1528 static void check_alloc_index(t_blocka *ba,int nindex_max)
1529 {
1530     if (nindex_max+1 > ba->nalloc_index)
1531     {
1532         ba->nalloc_index = over_alloc_dd(nindex_max+1);
1533         srenew(ba->index,ba->nalloc_index);
1534     }
1535 }
1536
1537 static void check_exclusions_alloc(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1538                                    t_blocka *lexcls)
1539 {
1540     int nr;
1541     int thread;
1542
1543     nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1544
1545     check_alloc_index(lexcls,nr);
1546
1547     for(thread=1; thread<dd->reverse_top->nthread; thread++)
1548     {
1549         check_alloc_index(&dd->reverse_top->excl_thread[thread],nr);
1550     }
1551 }
1552
1553 static void finish_local_exclusions(gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1554                                     t_blocka *lexcls)
1555 {
1556     int la0,la;
1557
1558     lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
1559
1560     if (dd->n_intercg_excl == 0)
1561     {
1562         /* There are no exclusions involving non-home charge groups,
1563          * but we need to set the indices for neighborsearching.
1564          */
1565         la0 = dd->cgindex[zones->izone[0].cg1];
1566         for(la=la0; la<lexcls->nr; la++)
1567         {
1568             lexcls->index[la] = lexcls->nra;
1569         }
1570
1571         /* nr is only used to loop over the exclusions for Ewald and RF,
1572          * so we can set it to the number of home atoms for efficiency.
1573          */
1574         lexcls->nr = dd->cgindex[zones->izone[0].cg1];
1575     }
1576 }
1577
1578 static void clear_idef(t_idef *idef)
1579 {
1580     int  ftype;
1581
1582      /* Clear the counts */
1583     for(ftype=0; ftype<F_NRE; ftype++)
1584     {
1585         idef->il[ftype].nr = 0;
1586     }
1587 }
1588
1589 static int make_local_bondeds_excls(gmx_domdec_t *dd,
1590                                     gmx_domdec_zones_t *zones,
1591                                     const gmx_mtop_t *mtop,
1592                                     const int *cginfo,
1593                                     gmx_bool bRCheckMB,ivec rcheck,gmx_bool bRCheck2B,
1594                                     real rc,
1595                                     int *la2lc,t_pbc *pbc_null,rvec *cg_cm,
1596                                     t_idef *idef,gmx_vsite_t *vsite,
1597                                     t_blocka *lexcls,int *excl_count)
1598 {
1599     int  nzone_bondeds,nzone_excl;
1600     int  iz,cg0,cg1;
1601     real rc2;
1602     int  nbonded_local;
1603     int  thread;
1604     gmx_reverse_top_t *rt;
1605
1606     if (dd->reverse_top->bMultiCGmols)
1607     {
1608         nzone_bondeds = zones->n;
1609     }
1610     else
1611     {
1612         /* Only single charge group molecules, so interactions don't
1613          * cross zone boundaries and we only need to assign in the home zone.
1614          */
1615         nzone_bondeds = 1;
1616     }
1617
1618     if (dd->n_intercg_excl > 0)
1619     {
1620         /* We only use exclusions from i-zones to i- and j-zones */
1621         nzone_excl = zones->nizone;
1622     }
1623     else
1624     {
1625         /* There are no inter-cg exclusions and only zone 0 sees itself */
1626         nzone_excl = 1;
1627     }
1628
1629     check_exclusions_alloc(dd,zones,lexcls);
1630
1631     rt = dd->reverse_top;
1632
1633     rc2 = rc*rc;
1634
1635     /* Clear the counts */
1636     clear_idef(idef);
1637     nbonded_local = 0;
1638
1639     lexcls->nr    = 0;
1640     lexcls->nra   = 0;
1641     *excl_count   = 0;
1642
1643     for(iz=0; iz<nzone_bondeds; iz++)
1644     {
1645         cg0 = zones->cg_range[iz];
1646         cg1 = zones->cg_range[iz+1];
1647
1648 #pragma omp parallel for num_threads(rt->nthread) schedule(static)
1649         for(thread=0; thread<rt->nthread; thread++)
1650         {
1651             int cg0t,cg1t;
1652             t_idef *idef_t;
1653             int ftype;
1654             int **vsite_pbc;
1655             int *vsite_pbc_nalloc;
1656             t_blocka *excl_t;
1657
1658             cg0t = cg0 + ((cg1 - cg0)* thread   )/rt->nthread;
1659             cg1t = cg0 + ((cg1 - cg0)*(thread+1))/rt->nthread;
1660
1661             if (thread == 0)
1662             {
1663                 idef_t = idef;
1664             }
1665             else
1666             {
1667                 idef_t = &rt->idef_thread[thread];
1668                 clear_idef(idef_t);
1669             }
1670
1671             if (vsite && vsite->bHaveChargeGroups && vsite->n_intercg_vsite > 0)
1672             {
1673                 if (thread == 0)
1674                 {
1675                     vsite_pbc        = vsite->vsite_pbc_loc;
1676                     vsite_pbc_nalloc = vsite->vsite_pbc_loc_nalloc;
1677                 }
1678                 else
1679                 {
1680                     vsite_pbc        = rt->vsite_pbc[thread];
1681                     vsite_pbc_nalloc = rt->vsite_pbc_nalloc[thread];
1682                 }
1683             }
1684             else
1685             {
1686                 vsite_pbc        = NULL;
1687                 vsite_pbc_nalloc = NULL;
1688             }
1689
1690             rt->nbonded_thread[thread] =
1691                 make_bondeds_zone(dd,zones,
1692                                   mtop->molblock,
1693                                   bRCheckMB,rcheck,bRCheck2B,rc2,
1694                                   la2lc,pbc_null,cg_cm,idef->iparams,
1695                                   idef_t,
1696                                   vsite,vsite_pbc,vsite_pbc_nalloc,
1697                                   iz,zones->n,
1698                                   dd->cgindex[cg0t],dd->cgindex[cg1t]);
1699
1700             if (iz < nzone_excl)
1701             {
1702                 if (thread == 0)
1703                 {
1704                     excl_t = lexcls;
1705                 }
1706                 else
1707                 {
1708                     excl_t = &rt->excl_thread[thread];
1709                     excl_t->nr  = 0;
1710                     excl_t->nra = 0;
1711                 }
1712
1713                 rt->excl_count_thread[thread] =
1714                     make_exclusions_zone(dd,zones,
1715                                          mtop->moltype,bRCheck2B,rc2,
1716                                          la2lc,pbc_null,cg_cm,cginfo,
1717                                          excl_t,
1718                                          iz,
1719                                          cg0t,cg1t);
1720             }
1721         }
1722
1723         if (rt->nthread > 1)
1724         {
1725             combine_idef(idef,rt->idef_thread+1,rt->nthread-1,
1726                          vsite,rt->vsite_pbc+1);
1727         }
1728
1729         for(thread=0; thread<rt->nthread; thread++)
1730         {
1731             nbonded_local += rt->nbonded_thread[thread];
1732         }
1733
1734         if (iz < nzone_excl)
1735         {
1736             if (rt->nthread > 1)
1737             {
1738                 combine_blocka(lexcls,rt->excl_thread+1,rt->nthread-1);
1739             }
1740
1741             for(thread=0; thread<rt->nthread; thread++)
1742             {
1743                 *excl_count += rt->excl_count_thread[thread];
1744             }
1745         }
1746     }
1747
1748     /* Some zones might not have exclusions, but some code still needs to
1749      * loop over the index, so we set the indices here.
1750      */
1751     for(iz=nzone_excl; iz<zones->nizone; iz++)
1752     {
1753         set_no_exclusions_zone(dd,zones,iz,lexcls);
1754     }
1755
1756     finish_local_exclusions(dd,zones,lexcls);
1757     if (debug)
1758     {
1759         fprintf(debug,"We have %d exclusions, check count %d\n",
1760                 lexcls->nra,*excl_count);
1761     }
1762
1763     return nbonded_local;
1764 }
1765
1766 void dd_make_local_cgs(gmx_domdec_t *dd,t_block *lcgs)
1767 {
1768   lcgs->nr    = dd->ncg_tot;
1769   lcgs->index = dd->cgindex;
1770 }
1771
1772 void dd_make_local_top(FILE *fplog,
1773                        gmx_domdec_t *dd,gmx_domdec_zones_t *zones,
1774                        int npbcdim,matrix box,
1775                        rvec cellsize_min,ivec npulse,
1776                        t_forcerec *fr,
1777                        rvec *cgcm_or_x,
1778                        gmx_vsite_t *vsite,
1779                        gmx_mtop_t *mtop,gmx_localtop_t *ltop)
1780 {
1781     gmx_bool bUniqueExcl,bRCheckMB,bRCheck2B,bRCheckExcl;
1782     real rc=-1;
1783     ivec rcheck;
1784     int  d,nexcl;
1785     t_pbc pbc,*pbc_null=NULL;
1786
1787     if (debug)
1788     {
1789         fprintf(debug,"Making local topology\n");
1790     }
1791
1792     dd_make_local_cgs(dd,&ltop->cgs);
1793
1794     bRCheckMB   = FALSE;
1795     bRCheck2B   = FALSE;
1796     bRCheckExcl = FALSE;
1797
1798     if (dd->reverse_top->bMultiCGmols)
1799     {
1800         /* We need to check to which cell bondeds should be assigned */
1801         rc = dd_cutoff_twobody(dd);
1802         if (debug)
1803         {
1804             fprintf(debug,"Two-body bonded cut-off distance is %g\n",rc);
1805         }
1806
1807         /* Should we check cg_cm distances when assigning bonded interactions? */
1808         for(d=0; d<DIM; d++)
1809         {
1810             rcheck[d] = FALSE;
1811             /* Only need to check for dimensions where the part of the box
1812              * that is not communicated is smaller than the cut-off.
1813              */
1814             if (d < npbcdim && dd->nc[d] > 1 &&
1815                 (dd->nc[d] - npulse[d])*cellsize_min[d] < 2*rc)
1816             {
1817                 if (dd->nc[d] == 2)
1818                 {
1819                     rcheck[d] = TRUE;
1820                     bRCheckMB = TRUE;
1821                 }
1822                 /* Check for interactions between two atoms,
1823                  * where we can allow interactions up to the cut-off,
1824                  * instead of up to the smallest cell dimension.
1825                  */
1826                 bRCheck2B = TRUE;
1827             }
1828             if (debug)
1829             {
1830                 fprintf(debug,
1831                         "dim %d cellmin %f bonded rcheck[%d] = %d, bRCheck2B = %d\n",
1832                         d,cellsize_min[d],d,rcheck[d],bRCheck2B);
1833             }
1834         }
1835         if (dd->reverse_top->bExclRequired)
1836         {
1837             bRCheckExcl = bRCheck2B;
1838         }
1839         else
1840         {
1841             /* If we don't have forces on exclusions,
1842              * we don't care about exclusions being assigned mulitple times.
1843              */
1844             bRCheckExcl = FALSE;
1845         }
1846         if (bRCheckMB || bRCheck2B)
1847         {
1848             make_la2lc(dd);
1849             if (fr->bMolPBC)
1850             {
1851                 set_pbc_dd(&pbc,fr->ePBC,dd,TRUE,box);
1852                 pbc_null = &pbc;
1853             }
1854             else
1855             {
1856                 pbc_null = NULL;
1857             }
1858         }
1859     }
1860
1861     dd->nbonded_local =
1862         make_local_bondeds_excls(dd,zones,mtop,fr->cginfo,
1863                                  bRCheckMB,rcheck,bRCheck2B,rc,
1864                                  dd->la2lc,
1865                                  pbc_null,cgcm_or_x,
1866                                  &ltop->idef,vsite,
1867                                  &ltop->excls,&nexcl);
1868
1869     /* The ilist is not sorted yet,
1870      * we can only do this when we have the charge arrays.
1871      */
1872     ltop->idef.ilsort = ilsortUNKNOWN;
1873
1874     if (dd->reverse_top->bExclRequired)
1875     {
1876         dd->nbonded_local += nexcl;
1877
1878         forcerec_set_excl_load(fr,ltop,NULL);
1879     }
1880
1881     ltop->atomtypes  = mtop->atomtypes;
1882
1883     /* For an error message only */
1884     dd->reverse_top->err_top_global = mtop;
1885     dd->reverse_top->err_top_local  = ltop;
1886 }
1887
1888 void dd_sort_local_top(gmx_domdec_t *dd,t_mdatoms *mdatoms,
1889                        gmx_localtop_t *ltop)
1890 {
1891     if (dd->reverse_top->ilsort == ilsortNO_FE)
1892     {
1893         ltop->idef.ilsort = ilsortNO_FE;
1894     }
1895     else
1896     {
1897         gmx_sort_ilist_fe(&ltop->idef,mdatoms->chargeA,mdatoms->chargeB);
1898     }
1899 }
1900
1901 gmx_localtop_t *dd_init_local_top(gmx_mtop_t *top_global)
1902 {
1903     gmx_localtop_t *top;
1904     int i;
1905
1906     snew(top,1);
1907
1908     top->idef.ntypes   = top_global->ffparams.ntypes;
1909     top->idef.atnr     = top_global->ffparams.atnr;
1910     top->idef.functype = top_global->ffparams.functype;
1911     top->idef.iparams  = top_global->ffparams.iparams;
1912     top->idef.fudgeQQ  = top_global->ffparams.fudgeQQ;
1913     top->idef.cmap_grid= top_global->ffparams.cmap_grid;
1914
1915     for(i=0; i<F_NRE; i++)
1916     {
1917         top->idef.il[i].iatoms = NULL;
1918         top->idef.il[i].nalloc = 0;
1919     }
1920     top->idef.ilsort   = ilsortUNKNOWN;
1921
1922     return top;
1923 }
1924
1925 void dd_init_local_state(gmx_domdec_t *dd,
1926                          t_state *state_global,t_state *state_local)
1927 {
1928     int buf[NITEM_DD_INIT_LOCAL_STATE];
1929
1930     if (DDMASTER(dd))
1931     {
1932         buf[0] = state_global->flags;
1933         buf[1] = state_global->ngtc;
1934         buf[2] = state_global->nnhpres;
1935         buf[3] = state_global->nhchainlength;
1936         buf[4] = state_global->dfhist.nlambda;
1937     }
1938     dd_bcast(dd,NITEM_DD_INIT_LOCAL_STATE*sizeof(int),buf);
1939
1940     init_state(state_local,0,buf[1],buf[2],buf[3],buf[4]);
1941     state_local->flags = buf[0];
1942
1943     /* With Langevin Dynamics we need to make proper storage space
1944      * in the global and local state for the random numbers.
1945      */
1946     if (state_local->flags & (1<<estLD_RNG))
1947     {
1948         if (DDMASTER(dd) && state_global->nrngi > 1)
1949         {
1950             state_global->nrng = dd->nnodes*gmx_rng_n();
1951             srenew(state_global->ld_rng,state_global->nrng);
1952         }
1953         state_local->nrng = gmx_rng_n();
1954         snew(state_local->ld_rng,state_local->nrng);
1955     }
1956     if (state_local->flags & (1<<estLD_RNGI))
1957     {
1958         if (DDMASTER(dd) && state_global->nrngi > 1)
1959         {
1960             state_global->nrngi = dd->nnodes;
1961             srenew(state_global->ld_rngi,state_global->nrngi);
1962         }
1963         snew(state_local->ld_rngi,1);
1964     }
1965 }
1966
1967 static void check_link(t_blocka *link,int cg_gl,int cg_gl_j)
1968 {
1969     int  k,aj;
1970     gmx_bool bFound;
1971
1972     bFound = FALSE;
1973     for(k=link->index[cg_gl]; k<link->index[cg_gl+1]; k++)
1974     {
1975         if (link->a[k] == cg_gl_j)
1976         {
1977             bFound = TRUE;
1978         }
1979     }
1980     if (!bFound)
1981     {
1982         /* Add this charge group link */
1983         if (link->index[cg_gl+1]+1 > link->nalloc_a)
1984         {
1985             link->nalloc_a = over_alloc_large(link->index[cg_gl+1]+1);
1986             srenew(link->a,link->nalloc_a);
1987         }
1988         link->a[link->index[cg_gl+1]] = cg_gl_j;
1989         link->index[cg_gl+1]++;
1990     }
1991 }
1992
1993 static int *make_at2cg(t_block *cgs)
1994 {
1995     int *at2cg,cg,a;
1996
1997     snew(at2cg,cgs->index[cgs->nr]);
1998     for(cg=0; cg<cgs->nr; cg++)
1999     {
2000         for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
2001         {
2002             at2cg[a] = cg;
2003         }
2004     }
2005
2006     return at2cg;
2007 }
2008
2009 t_blocka *make_charge_group_links(gmx_mtop_t *mtop,gmx_domdec_t *dd,
2010                                   cginfo_mb_t *cginfo_mb)
2011 {
2012     gmx_reverse_top_t *rt;
2013     int  mb,cg_offset,cg,cg_gl,a,aj,i,j,ftype,nral,nlink_mol,mol,ncgi;
2014     gmx_molblock_t *molb;
2015     gmx_moltype_t *molt;
2016     t_block *cgs;
2017     t_blocka *excls;
2018     int *a2c;
2019     gmx_reverse_ilist_t ril;
2020     t_blocka *link;
2021     cginfo_mb_t *cgi_mb;
2022
2023     /* For each charge group make a list of other charge groups
2024      * in the system that a linked to it via bonded interactions
2025      * which are also stored in reverse_top.
2026      */
2027
2028     rt = dd->reverse_top;
2029
2030     snew(link,1);
2031     snew(link->index,ncg_mtop(mtop)+1);
2032     link->nalloc_a = 0;
2033     link->a = NULL;
2034
2035     link->index[0] = 0;
2036     cg_offset = 0;
2037     ncgi = 0;
2038     for(mb=0; mb<mtop->nmolblock; mb++)
2039     {
2040         molb = &mtop->molblock[mb];
2041         if (molb->nmol == 0)
2042         {
2043             continue;
2044         }
2045         molt = &mtop->moltype[molb->type];
2046         cgs   = &molt->cgs;
2047         excls = &molt->excls;
2048         a2c = make_at2cg(cgs);
2049         /* Make a reverse ilist in which the interactions are linked
2050          * to all atoms, not only the first atom as in gmx_reverse_top.
2051          * The constraints are discarded here.
2052          */
2053         make_reverse_ilist(molt,NULL,FALSE,FALSE,FALSE,TRUE,&ril);
2054
2055         cgi_mb = &cginfo_mb[mb];
2056
2057         for(cg=0; cg<cgs->nr; cg++)
2058         {
2059             cg_gl = cg_offset + cg;
2060             link->index[cg_gl+1] = link->index[cg_gl];
2061             for(a=cgs->index[cg]; a<cgs->index[cg+1]; a++)
2062             {
2063                 i = ril.index[a];
2064                 while (i < ril.index[a+1])
2065                 {
2066                     ftype = ril.il[i++];
2067                     nral = NRAL(ftype);
2068                     /* Skip the ifunc index */
2069                     i++;
2070                     for(j=0; j<nral; j++)
2071                     {
2072                         aj = ril.il[i+j];
2073                         if (a2c[aj] != cg)
2074                         {
2075                             check_link(link,cg_gl,cg_offset+a2c[aj]);
2076                         }
2077                     }
2078                     i += nral_rt(ftype);
2079                 }
2080                 if (rt->bExclRequired)
2081                 {
2082                     /* Exclusions always go both ways */
2083                     for(j=excls->index[a]; j<excls->index[a+1]; j++)
2084                     {
2085                         aj = excls->a[j];
2086                         if (a2c[aj] != cg)
2087                         {
2088                             check_link(link,cg_gl,cg_offset+a2c[aj]);
2089                         }
2090                     }
2091                 }
2092             }
2093             if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
2094             {
2095                 SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
2096                 ncgi++;
2097             }
2098         }
2099         nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
2100
2101         cg_offset += cgs->nr;
2102
2103         destroy_reverse_ilist(&ril);
2104         sfree(a2c);
2105
2106         if (debug)
2107         {
2108             fprintf(debug,"molecule type '%s' %d cgs has %d cg links through bonded interac.\n",*molt->name,cgs->nr,nlink_mol);
2109         }
2110
2111         if (molb->nmol > 1)
2112         {
2113             /* Copy the data for the rest of the molecules in this block */
2114             link->nalloc_a += (molb->nmol - 1)*nlink_mol;
2115             srenew(link->a,link->nalloc_a);
2116             for(mol=1; mol<molb->nmol; mol++)
2117             {
2118                 for(cg=0; cg<cgs->nr; cg++)
2119                 {
2120                     cg_gl = cg_offset + cg;
2121                     link->index[cg_gl+1] =
2122                         link->index[cg_gl+1-cgs->nr] + nlink_mol;
2123                     for(j=link->index[cg_gl]; j<link->index[cg_gl+1]; j++)
2124                     {
2125                         link->a[j] = link->a[j-nlink_mol] + cgs->nr;
2126                     }
2127                     if (link->index[cg_gl+1] - link->index[cg_gl] > 0 &&
2128                         cg_gl - cgi_mb->cg_start < cgi_mb->cg_mod)
2129                     {
2130                         SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg_gl - cgi_mb->cg_start]);
2131                         ncgi++;
2132                     }
2133                 }
2134                 cg_offset += cgs->nr;
2135             }
2136         }
2137     }
2138
2139     if (debug)
2140     {
2141         fprintf(debug,"Of the %d charge groups %d are linked via bonded interactions\n",ncg_mtop(mtop),ncgi);
2142     }
2143
2144     return link;
2145 }
2146
2147 static void bonded_cg_distance_mol(gmx_moltype_t *molt,int *at2cg,
2148                                    gmx_bool bBCheck,gmx_bool bExcl,rvec *cg_cm,
2149                                    real *r_2b,int *ft2b,int *a2_1,int *a2_2,
2150                                    real *r_mb,int *ftmb,int *am_1,int *am_2)
2151 {
2152     int ftype,nral,i,j,ai,aj,cgi,cgj;
2153     t_ilist *il;
2154     t_blocka *excls;
2155     real r2_2b,r2_mb,rij2;
2156
2157     r2_2b = 0;
2158     r2_mb = 0;
2159     for(ftype=0; ftype<F_NRE; ftype++)
2160     {
2161         if (dd_check_ftype(ftype,bBCheck,FALSE,FALSE))
2162         {
2163             il = &molt->ilist[ftype];
2164             nral = NRAL(ftype);
2165             if (nral > 1)
2166             {
2167                 for(i=0; i<il->nr; i+=1+nral)
2168                 {
2169                     for(ai=0; ai<nral; ai++)
2170                     {
2171                         cgi = at2cg[il->iatoms[i+1+ai]];
2172                                                 for(aj=0; aj<nral; aj++) {
2173                                                         cgj = at2cg[il->iatoms[i+1+aj]];
2174                                                         if (cgi != cgj)
2175                             {
2176                                                                 rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
2177                                                                 if (nral == 2 && rij2 > r2_2b)
2178                                 {
2179                                                                         r2_2b = rij2;
2180                                     *ft2b = ftype;
2181                                     *a2_1 = il->iatoms[i+1+ai];
2182                                     *a2_2 = il->iatoms[i+1+aj];
2183                                                                 }
2184                                                                 if (nral >  2 && rij2 > r2_mb)
2185                                 {
2186                                                                         r2_mb = rij2;
2187                                     *ftmb = ftype;
2188                                     *am_1 = il->iatoms[i+1+ai];
2189                                     *am_2 = il->iatoms[i+1+aj];
2190                                                                 }
2191                                                         }
2192                                                 }
2193                                         }
2194                                 }
2195                         }
2196                 }
2197         }
2198         if (bExcl)
2199     {
2200                 excls = &molt->excls;
2201                 for(ai=0; ai<excls->nr; ai++)
2202         {
2203                         cgi = at2cg[ai];
2204                         for(j=excls->index[ai]; j<excls->index[ai+1]; j++) {
2205                                 cgj = at2cg[excls->a[j]];
2206                                 if (cgi != cgj)
2207                 {
2208                                         rij2 = distance2(cg_cm[cgi],cg_cm[cgj]);
2209                                         if (rij2 > r2_2b)
2210                     {
2211                                                 r2_2b = rij2;
2212                                         }
2213                                 }
2214                         }
2215                 }
2216         }
2217
2218         *r_2b = sqrt(r2_2b);
2219         *r_mb = sqrt(r2_mb);
2220 }
2221
2222 static void get_cgcm_mol(gmx_moltype_t *molt,gmx_ffparams_t *ffparams,
2223                          int ePBC,t_graph *graph,matrix box,
2224                          gmx_vsite_t *vsite,
2225                          rvec *x,rvec *xs,rvec *cg_cm)
2226 {
2227     int n,i;
2228
2229     if (ePBC != epbcNONE)
2230     {
2231         mk_mshift(NULL,graph,ePBC,box,x);
2232
2233         shift_x(graph,box,x,xs);
2234         /* By doing an extra mk_mshift the molecules that are broken
2235          * because they were e.g. imported from another software
2236          * will be made whole again. Such are the healing powers
2237          * of GROMACS.
2238          */
2239         mk_mshift(NULL,graph,ePBC,box,xs);
2240     }
2241     else
2242     {
2243         /* We copy the coordinates so the original coordinates remain
2244          * unchanged, just to be 100% sure that we do not affect
2245          * binary reproducibility of simulations.
2246          */
2247         n = molt->cgs.index[molt->cgs.nr];
2248         for(i=0; i<n; i++)
2249         {
2250             copy_rvec(x[i],xs[i]);
2251         }
2252     }
2253
2254     if (vsite)
2255     {
2256         construct_vsites(NULL,vsite,xs,NULL,0.0,NULL,
2257                          ffparams->iparams,molt->ilist,
2258                          epbcNONE,TRUE,NULL,NULL,NULL);
2259     }
2260
2261     calc_cgcm(NULL,0,molt->cgs.nr,&molt->cgs,xs,cg_cm);
2262 }
2263
2264 static int have_vsite_molt(gmx_moltype_t *molt)
2265 {
2266     int  i;
2267     gmx_bool bVSite;
2268
2269     bVSite = FALSE;
2270     for(i=0; i<F_NRE; i++)
2271     {
2272         if ((interaction_function[i].flags & IF_VSITE) &&
2273             molt->ilist[i].nr > 0) {
2274             bVSite = TRUE;
2275         }
2276     }
2277
2278     return bVSite;
2279 }
2280
2281 void dd_bonded_cg_distance(FILE *fplog,
2282                            gmx_domdec_t *dd,gmx_mtop_t *mtop,
2283                            t_inputrec *ir,rvec *x,matrix box,
2284                            gmx_bool bBCheck,
2285                            real *r_2b,real *r_mb)
2286 {
2287     gmx_bool bExclRequired;
2288     int  mb,cg_offset,at_offset,*at2cg,mol;
2289     t_graph graph;
2290     gmx_vsite_t *vsite;
2291     gmx_molblock_t *molb;
2292     gmx_moltype_t *molt;
2293     rvec *xs,*cg_cm;
2294     real rmol_2b,rmol_mb;
2295     int ft2b=-1,a_2b_1=-1,a_2b_2=-1,ftmb=-1,a_mb_1=-1,a_mb_2=-1;
2296     int ftm2b=-1,amol_2b_1=-1,amol_2b_2=-1,ftmmb=-1,amol_mb_1=-1,amol_mb_2=-1;
2297
2298     bExclRequired = IR_EXCL_FORCES(*ir);
2299
2300     vsite = init_vsite(mtop,NULL,TRUE);
2301
2302     *r_2b = 0;
2303     *r_mb = 0;
2304     cg_offset = 0;
2305     at_offset = 0;
2306     for(mb=0; mb<mtop->nmolblock; mb++)
2307     {
2308         molb = &mtop->molblock[mb];
2309         molt = &mtop->moltype[molb->type];
2310         if (molt->cgs.nr == 1 || molb->nmol == 0)
2311         {
2312             cg_offset += molb->nmol*molt->cgs.nr;
2313             at_offset += molb->nmol*molt->atoms.nr;
2314         }
2315         else
2316         {
2317             if (ir->ePBC != epbcNONE)
2318             {
2319                 mk_graph_ilist(NULL,molt->ilist,0,molt->atoms.nr,FALSE,FALSE,
2320                                &graph);
2321             }
2322
2323             at2cg = make_at2cg(&molt->cgs);
2324             snew(xs,molt->atoms.nr);
2325             snew(cg_cm,molt->cgs.nr);
2326             for(mol=0; mol<molb->nmol; mol++)
2327             {
2328                 get_cgcm_mol(molt,&mtop->ffparams,ir->ePBC,&graph,box,
2329                              have_vsite_molt(molt) ? vsite : NULL,
2330                              x+at_offset,xs,cg_cm);
2331
2332                 bonded_cg_distance_mol(molt,at2cg,bBCheck,bExclRequired,cg_cm,
2333                                        &rmol_2b,&ftm2b,&amol_2b_1,&amol_2b_2,
2334                                        &rmol_mb,&ftmmb,&amol_mb_1,&amol_mb_2);
2335                 if (rmol_2b > *r_2b)
2336                 {
2337                     *r_2b  = rmol_2b;
2338                     ft2b   = ftm2b;
2339                     a_2b_1 = at_offset + amol_2b_1;
2340                     a_2b_2 = at_offset + amol_2b_2;
2341                 }
2342                 if (rmol_mb > *r_mb)
2343                 {
2344                     *r_mb  = rmol_mb;
2345                     ftmb   = ftmmb;
2346                     a_mb_1 = at_offset + amol_mb_1;
2347                     a_mb_2 = at_offset + amol_mb_2;
2348                 }
2349
2350                 cg_offset += molt->cgs.nr;
2351                 at_offset += molt->atoms.nr;
2352             }
2353             sfree(cg_cm);
2354             sfree(xs);
2355             sfree(at2cg);
2356             if (ir->ePBC != epbcNONE)
2357             {
2358                 done_graph(&graph);
2359             }
2360         }
2361     }
2362
2363     /* We should have a vsite free routine, but here we can simply free */
2364     sfree(vsite);
2365
2366     if (fplog && (ft2b >= 0 || ftmb >= 0))
2367     {
2368         fprintf(fplog,
2369                 "Initial maximum inter charge-group distances:\n");
2370         if (ft2b >= 0)
2371         {
2372             fprintf(fplog,
2373                     "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2374                     *r_2b,interaction_function[ft2b].longname,
2375                     a_2b_1+1,a_2b_2+1);
2376         }
2377         if (ftmb >= 0)
2378         {
2379             fprintf(fplog,
2380                     "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d\n",
2381                     *r_mb,interaction_function[ftmb].longname,
2382                     a_mb_1+1,a_mb_2+1);
2383         }
2384     }
2385 }