src/gromacs/mdlib/domdec_con.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2006,2007,2008,2009,2010,2012,2013, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39 #include <assert.h>
  40
  41 #include "smalloc.h"
  42 #include "vec.h"
  43 #include "constr.h"
  44 #include "domdec.h"
  45 #include "domdec_network.h"
  46 #include "mtop_util.h"
  47 #include "gmx_ga2la.h"
  48 #include "gmx_hash.h"
  49 #include "gmx_omp_nthreads.h"
  50 #include "macros.h"
  51
  52 typedef struct {
  53     int  nsend;
  54     int *a;
  55     int  a_nalloc;
  56     int  nrecv;
  57 } gmx_specatsend_t;
  58
  59 typedef struct {
  60     int *ind;
  61     int  nalloc;
  62     int  n;
  63 } ind_req_t;
  64
  65 typedef struct gmx_domdec_specat_comm {
  66     /* The number of indices to receive during the setup */
  67     int              nreq[DIM][2][2];
  68     /* The atoms to send */
  69     gmx_specatsend_t spas[DIM][2];
  70     gmx_bool        *bSendAtom;
  71     int              bSendAtom_nalloc;
  72     /* Send buffers */
  73     int             *ibuf;
  74     int              ibuf_nalloc;
  75     rvec            *vbuf;
  76     int              vbuf_nalloc;
  77     rvec            *vbuf2;
  78     int              vbuf2_nalloc;
  79     /* The range in the local buffer(s) for received atoms */
  80     int              at_start;
  81     int              at_end;
  82
  83     /* The atom indices we need from the surrounding cells.
  84      * We can gather the indices over nthread threads.
  85      */
  86     int        nthread;
  87     ind_req_t *ireq;
  88 } gmx_domdec_specat_comm_t;
  89
  90 typedef struct gmx_domdec_constraints {
  91     int       *molb_con_offset;
  92     int       *molb_ncon_mol;
  93     /* The fully local and connected constraints */
  94     int        ncon;
  95     /* The global constraint number, only required for clearing gc_req */
  96     int       *con_gl;
  97     int       *con_nlocat;
  98     int        con_nalloc;
  99     /* Boolean that tells if a global constraint index has been requested */
 100     char      *gc_req;
 101     /* Global to local communicated constraint atom only index */
 102     gmx_hash_t ga2la;
 103
 104     /* Multi-threading stuff */
 105     int      nthread;
 106     t_ilist *ils;
 107 } gmx_domdec_constraints_t;
 108
 109
 110 static void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
 111                              rvec *f, rvec *fshift)
 112 {
 113     gmx_specatsend_t *spas;
 114     rvec             *vbuf;
 115     int               n, n0, n1, d, dim, dir, i;
 116     ivec              vis;
 117     int               is;
 118     gmx_bool          bPBC, bScrew;
 119
 120     n = spac->at_end;
 121     for (d = dd->ndim-1; d >= 0; d--)
 122     {
 123         dim = dd->dim[d];
 124         if (dd->nc[dim] > 2)
 125         {
 126             /* Pulse the grid forward and backward */
 127             spas = spac->spas[d];
 128             n0   = spas[0].nrecv;
 129             n1   = spas[1].nrecv;
 130             n   -= n1 + n0;
 131             vbuf = spac->vbuf;
 132             /* Send and receive the coordinates */
 133             dd_sendrecv2_rvec(dd, d,
 134                               f+n+n1, n0, vbuf, spas[0].nsend,
 135                               f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
 136             for (dir = 0; dir < 2; dir++)
 137             {
 138                 bPBC   = ((dir == 0 && dd->ci[dim] == 0) ||
 139                           (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
 140                 bScrew = (bPBC && dd->bScrewPBC && dim == XX);
 141
 142                 spas = &spac->spas[d][dir];
 143                 /* Sum the buffer into the required forces */
 144                 if (!bPBC || (!bScrew && fshift == NULL))
 145                 {
 146                     for (i = 0; i < spas->nsend; i++)
 147                     {
 148                         rvec_inc(f[spas->a[i]], *vbuf);
 149                         vbuf++;
 150                     }
 151                 }
 152                 else
 153                 {
 154                     clear_ivec(vis);
 155                     vis[dim] = (dir == 0 ? 1 : -1);
 156                     is       = IVEC2IS(vis);
 157                     if (!bScrew)
 158                     {
 159                         /* Sum and add to shift forces */
 160                         for (i = 0; i < spas->nsend; i++)
 161                         {
 162                             rvec_inc(f[spas->a[i]], *vbuf);
 163                             rvec_inc(fshift[is], *vbuf);
 164                             vbuf++;
 165                         }
 166                     }
 167                     else
 168                     {
 169                         /* Rotate the forces */
 170                         for (i = 0; i < spas->nsend; i++)
 171                         {
 172                             f[spas->a[i]][XX] += (*vbuf)[XX];
 173                             f[spas->a[i]][YY] -= (*vbuf)[YY];
 174                             f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
 175                             if (fshift)
 176                             {
 177                                 rvec_inc(fshift[is], *vbuf);
 178                             }
 179                             vbuf++;
 180                         }
 181                     }
 182                 }
 183             }
 184         }
 185         else
 186         {
 187             /* Two cells, so we only need to communicate one way */
 188             spas = &spac->spas[d][0];
 189             n   -= spas->nrecv;
 190             /* Send and receive the coordinates */
 191             dd_sendrecv_rvec(dd, d, dddirForward,
 192                              f+n, spas->nrecv, spac->vbuf, spas->nsend);
 193             /* Sum the buffer into the required forces */
 194             if (dd->bScrewPBC && dim == XX &&
 195                 (dd->ci[dim] == 0 ||
 196                  dd->ci[dim] == dd->nc[dim]-1))
 197             {
 198                 for (i = 0; i < spas->nsend; i++)
 199                 {
 200                     /* Rotate the force */
 201                     f[spas->a[i]][XX] += spac->vbuf[i][XX];
 202                     f[spas->a[i]][YY] -= spac->vbuf[i][YY];
 203                     f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
 204                 }
 205             }
 206             else
 207             {
 208                 for (i = 0; i < spas->nsend; i++)
 209                 {
 210                     rvec_inc(f[spas->a[i]], spac->vbuf[i]);
 211                 }
 212             }
 213         }
 214     }
 215 }
 216
 217 void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
 218 {
 219     if (dd->vsite_comm)
 220     {
 221         dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
 222     }
 223 }
 224
 225 void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
 226 {
 227     int i;
 228
 229     if (dd->vsite_comm)
 230     {
 231         for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
 232         {
 233             clear_rvec(f[i]);
 234         }
 235     }
 236 }
 237
 238 static void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
 239                              matrix box, rvec *x0, rvec *x1)
 240 {
 241     gmx_specatsend_t *spas;
 242     rvec             *x, *vbuf, *rbuf;
 243     int               nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
 244     gmx_bool          bPBC, bScrew = FALSE;
 245     rvec              shift = {0, 0, 0};
 246
 247     nvec = 1;
 248     if (x1)
 249     {
 250         nvec++;
 251     }
 252
 253     n = spac->at_start;
 254     for (d = 0; d < dd->ndim; d++)
 255     {
 256         dim = dd->dim[d];
 257         if (dd->nc[dim] > 2)
 258         {
 259             /* Pulse the grid forward and backward */
 260             vbuf = spac->vbuf;
 261             for (dir = 0; dir < 2; dir++)
 262             {
 263                 if (dir == 0 && dd->ci[dim] == 0)
 264                 {
 265                     bPBC   = TRUE;
 266                     bScrew = (dd->bScrewPBC && dim == XX);
 267                     copy_rvec(box[dim], shift);
 268                 }
 269                 else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
 270                 {
 271                     bPBC   = TRUE;
 272                     bScrew = (dd->bScrewPBC && dim == XX);
 273                     for (i = 0; i < DIM; i++)
 274                     {
 275                         shift[i] = -box[dim][i];
 276                     }
 277                 }
 278                 else
 279                 {
 280                     bPBC   = FALSE;
 281                     bScrew = FALSE;
 282                 }
 283                 spas = &spac->spas[d][dir];
 284                 for (v = 0; v < nvec; v++)
 285                 {
 286                     x = (v == 0 ? x0 : x1);
 287                     /* Copy the required coordinates to the send buffer */
 288                     if (!bPBC)
 289                     {
 290                         /* Only copy */
 291                         for (i = 0; i < spas->nsend; i++)
 292                         {
 293                             copy_rvec(x[spas->a[i]], *vbuf);
 294                             vbuf++;
 295                         }
 296                     }
 297                     else if (!bScrew)
 298                     {
 299                         /* Shift coordinates */
 300                         for (i = 0; i < spas->nsend; i++)
 301                         {
 302                             rvec_add(x[spas->a[i]], shift, *vbuf);
 303                             vbuf++;
 304                         }
 305                     }
 306                     else
 307                     {
 308                         /* Shift and rotate coordinates */
 309                         for (i = 0; i < spas->nsend; i++)
 310                         {
 311                             (*vbuf)[XX] =               x[spas->a[i]][XX] + shift[XX];
 312                             (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
 313                             (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
 314                             vbuf++;
 315                         }
 316                     }
 317                 }
 318             }
 319             /* Send and receive the coordinates */
 320             spas = spac->spas[d];
 321             ns0  = spas[0].nsend;
 322             nr0  = spas[0].nrecv;
 323             ns1  = spas[1].nsend;
 324             nr1  = spas[1].nrecv;
 325             if (nvec == 1)
 326             {
 327                 dd_sendrecv2_rvec(dd, d,
 328                                   spac->vbuf+ns0, ns1, x0+n, nr1,
 329                                   spac->vbuf, ns0, x0+n+nr1, nr0);
 330             }
 331             else
 332             {
 333                 /* Communicate both vectors in one buffer */
 334                 rbuf = spac->vbuf2;
 335                 dd_sendrecv2_rvec(dd, d,
 336                                   spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
 337                                   spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
 338                 /* Split the buffer into the two vectors */
 339                 nn = n;
 340                 for (dir = 1; dir >= 0; dir--)
 341                 {
 342                     nr = spas[dir].nrecv;
 343                     for (v = 0; v < 2; v++)
 344                     {
 345                         x = (v == 0 ? x0 : x1);
 346                         for (i = 0; i < nr; i++)
 347                         {
 348                             copy_rvec(*rbuf, x[nn+i]);
 349                             rbuf++;
 350                         }
 351                     }
 352                     nn += nr;
 353                 }
 354             }
 355             n += nr0 + nr1;
 356         }
 357         else
 358         {
 359             spas = &spac->spas[d][0];
 360             /* Copy the required coordinates to the send buffer */
 361             vbuf = spac->vbuf;
 362             for (v = 0; v < nvec; v++)
 363             {
 364                 x = (v == 0 ? x0 : x1);
 365                 if (dd->bScrewPBC && dim == XX &&
 366                     (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
 367                 {
 368                     /* Here we only perform the rotation, the rest of the pbc
 369                      * is handled in the constraint or viste routines.
 370                      */
 371                     for (i = 0; i < spas->nsend; i++)
 372                     {
 373                         (*vbuf)[XX] =               x[spas->a[i]][XX];
 374                         (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
 375                         (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
 376                         vbuf++;
 377                     }
 378                 }
 379                 else
 380                 {
 381                     for (i = 0; i < spas->nsend; i++)
 382                     {
 383                         copy_rvec(x[spas->a[i]], *vbuf);
 384                         vbuf++;
 385                     }
 386                 }
 387             }
 388             /* Send and receive the coordinates */
 389             if (nvec == 1)
 390             {
 391                 dd_sendrecv_rvec(dd, d, dddirBackward,
 392                                  spac->vbuf, spas->nsend, x0+n, spas->nrecv);
 393             }
 394             else
 395             {
 396                 /* Communicate both vectors in one buffer */
 397                 rbuf = spac->vbuf2;
 398                 dd_sendrecv_rvec(dd, d, dddirBackward,
 399                                  spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
 400                 /* Split the buffer into the two vectors */
 401                 nr = spas[0].nrecv;
 402                 for (v = 0; v < 2; v++)
 403                 {
 404                     x = (v == 0 ? x0 : x1);
 405                     for (i = 0; i < nr; i++)
 406                     {
 407                         copy_rvec(*rbuf, x[n+i]);
 408                         rbuf++;
 409                     }
 410                 }
 411             }
 412             n += spas->nrecv;
 413         }
 414     }
 415 }
 416
 417 void dd_move_x_constraints(gmx_domdec_t *dd, matrix box, rvec *x0, rvec *x1)
 418 {
 419     if (dd->constraint_comm)
 420     {
 421         dd_move_x_specat(dd, dd->constraint_comm, box, x0, x1);
 422     }
 423 }
 424
 425 void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x)
 426 {
 427     if (dd->vsite_comm)
 428     {
 429         dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL);
 430     }
 431 }
 432
 433 int *dd_constraints_nlocalatoms(gmx_domdec_t *dd)
 434 {
 435     if (dd->constraints)
 436     {
 437         return dd->constraints->con_nlocat;
 438     }
 439     else
 440     {
 441         return NULL;
 442     }
 443 }
 444
 445 void dd_clear_local_constraint_indices(gmx_domdec_t *dd)
 446 {
 447     gmx_domdec_constraints_t *dc;
 448     int i;
 449
 450     dc = dd->constraints;
 451
 452     for (i = 0; i < dc->ncon; i++)
 453     {
 454         dc->gc_req[dc->con_gl[i]] = 0;
 455     }
 456
 457     if (dd->constraint_comm)
 458     {
 459         gmx_hash_clear_and_optimize(dc->ga2la);
 460     }
 461 }
 462
 463 void dd_clear_local_vsite_indices(gmx_domdec_t *dd)
 464 {
 465     int i;
 466
 467     if (dd->vsite_comm)
 468     {
 469         gmx_hash_clear_and_optimize(dd->ga2la_vsite);
 470     }
 471 }
 472
 473 static int setup_specat_communication(gmx_domdec_t             *dd,
 474                                       ind_req_t                *ireq,
 475                                       gmx_domdec_specat_comm_t *spac,
 476                                       gmx_hash_t                ga2la_specat,
 477                                       int                       at_start,
 478                                       int                       vbuf_fac,
 479                                       const char               *specat_type,
 480                                       const char               *add_err)
 481 {
 482     int               nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
 483     int               d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
 484     int               nat_tot_specat, nat_tot_prev, nalloc_old;
 485     gmx_bool          bPBC, bFirst;
 486     gmx_specatsend_t *spas;
 487
 488     if (debug)
 489     {
 490         fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
 491     }
 492
 493     /* nsend[0]: the number of atoms requested by this node only,
 494      *           we communicate this for more efficients checks
 495      * nsend[1]: the total number of requested atoms
 496      */
 497     nsend[0] = ireq->n;
 498     nsend[1] = nsend[0];
 499     nlast    = nsend[1];
 500     for (d = dd->ndim-1; d >= 0; d--)
 501     {
 502         /* Pulse the grid forward and backward */
 503         dim  = dd->dim[d];
 504         bPBC = (dim < dd->npbcdim);
 505         if (dd->nc[dim] == 2)
 506         {
 507             /* Only 2 cells, so we only need to communicate once */
 508             ndir = 1;
 509         }
 510         else
 511         {
 512             ndir = 2;
 513         }
 514         for (dir = 0; dir < ndir; dir++)
 515         {
 516             if (!bPBC &&
 517                 dd->nc[dim] > 2 &&
 518                 ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
 519                  (dir == 1 && dd->ci[dim] == 0)))
 520             {
 521                 /* No pbc: the fist/last cell should not request atoms */
 522                 nsend_ptr = nsend_zero;
 523             }
 524             else
 525             {
 526                 nsend_ptr = nsend;
 527             }
 528             /* Communicate the number of indices */
 529             dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
 530                             nsend_ptr, 2, spac->nreq[d][dir], 2);
 531             nr = spac->nreq[d][dir][1];
 532             if (nlast+nr > ireq->nalloc)
 533             {
 534                 ireq->nalloc = over_alloc_dd(nlast+nr);
 535                 srenew(ireq->ind, ireq->nalloc);
 536             }
 537             /* Communicate the indices */
 538             dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
 539                             ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
 540             nlast += nr;
 541         }
 542         nsend[1] = nlast;
 543     }
 544     if (debug)
 545     {
 546         fprintf(debug, "Communicated the counts\n");
 547     }
 548
 549     /* Search for the requested atoms and communicate the indices we have */
 550     nat_tot_specat = at_start;
 551     nrecv_local    = 0;
 552     for (d = 0; d < dd->ndim; d++)
 553     {
 554         bFirst = (d == 0);
 555         /* Pulse the grid forward and backward */
 556         if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
 557         {
 558             ndir = 2;
 559         }
 560         else
 561         {
 562             ndir = 1;
 563         }
 564         nat_tot_prev = nat_tot_specat;
 565         for (dir = ndir-1; dir >= 0; dir--)
 566         {
 567             if (nat_tot_specat > spac->bSendAtom_nalloc)
 568             {
 569                 nalloc_old             = spac->bSendAtom_nalloc;
 570                 spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
 571                 srenew(spac->bSendAtom, spac->bSendAtom_nalloc);
 572                 for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++)
 573                 {
 574                     spac->bSendAtom[i] = FALSE;
 575                 }
 576             }
 577             spas = &spac->spas[d][dir];
 578             n0   = spac->nreq[d][dir][0];
 579             nr   = spac->nreq[d][dir][1];
 580             if (debug)
 581             {
 582                 fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n",
 583                         d, dir, nr);
 584             }
 585             start       = nlast - nr;
 586             spas->nsend = 0;
 587             nsend[0]    = 0;
 588             for (i = 0; i < nr; i++)
 589             {
 590                 indr = ireq->ind[start+i];
 591                 ind  = -1;
 592                 /* Check if this is a home atom and if so ind will be set */
 593                 if (!ga2la_get_home(dd->ga2la, indr, &ind))
 594                 {
 595                     /* Search in the communicated atoms */
 596                     ind = gmx_hash_get_minone(ga2la_specat, indr);
 597                 }
 598                 if (ind >= 0)
 599                 {
 600                     if (i < n0 || !spac->bSendAtom[ind])
 601                     {
 602                         if (spas->nsend+1 > spas->a_nalloc)
 603                         {
 604                             spas->a_nalloc = over_alloc_large(spas->nsend+1);
 605                             srenew(spas->a, spas->a_nalloc);
 606                         }
 607                         /* Store the local index so we know which coordinates
 608                          * to send out later.
 609                          */
 610                         spas->a[spas->nsend] = ind;
 611                         spac->bSendAtom[ind] = TRUE;
 612                         if (spas->nsend+1 > spac->ibuf_nalloc)
 613                         {
 614                             spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
 615                             srenew(spac->ibuf, spac->ibuf_nalloc);
 616                         }
 617                         /* Store the global index so we can send it now */
 618                         spac->ibuf[spas->nsend] = indr;
 619                         if (i < n0)
 620                         {
 621                             nsend[0]++;
 622                         }
 623                         spas->nsend++;
 624                     }
 625                 }
 626             }
 627             nlast = start;
 628             /* Clear the local flags */
 629             for (i = 0; i < spas->nsend; i++)
 630             {
 631                 spac->bSendAtom[spas->a[i]] = FALSE;
 632             }
 633             /* Send and receive the number of indices to communicate */
 634             nsend[1] = spas->nsend;
 635             dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
 636                             nsend, 2, buf, 2);
 637             if (debug)
 638             {
 639                 fprintf(debug, "Send to node %d, %d (%d) indices, "
 640                         "receive from node %d, %d (%d) indices\n",
 641                         dd->neighbor[d][1-dir], nsend[1], nsend[0],
 642                         dd->neighbor[d][dir], buf[1], buf[0]);
 643                 if (gmx_debug_at)
 644                 {
 645                     for (i = 0; i < spas->nsend; i++)
 646                     {
 647                         fprintf(debug, " %d", spac->ibuf[i]+1);
 648                     }
 649                     fprintf(debug, "\n");
 650                 }
 651             }
 652             nrecv_local += buf[0];
 653             spas->nrecv  = buf[1];
 654             if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
 655             {
 656                 dd->gatindex_nalloc =
 657                     over_alloc_dd(nat_tot_specat + spas->nrecv);
 658                 srenew(dd->gatindex, dd->gatindex_nalloc);
 659             }
 660             /* Send and receive the indices */
 661             dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
 662                             spac->ibuf, spas->nsend,
 663                             dd->gatindex+nat_tot_specat, spas->nrecv);
 664             nat_tot_specat += spas->nrecv;
 665         }
 666
 667         /* Allocate the x/f communication buffers */
 668         ns = spac->spas[d][0].nsend;
 669         nr = spac->spas[d][0].nrecv;
 670         if (ndir == 2)
 671         {
 672             ns += spac->spas[d][1].nsend;
 673             nr += spac->spas[d][1].nrecv;
 674         }
 675         if (vbuf_fac*ns > spac->vbuf_nalloc)
 676         {
 677             spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
 678             srenew(spac->vbuf, spac->vbuf_nalloc);
 679         }
 680         if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
 681         {
 682             spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
 683             srenew(spac->vbuf2, spac->vbuf2_nalloc);
 684         }
 685
 686         /* Make a global to local index for the communication atoms */
 687         for (i = nat_tot_prev; i < nat_tot_specat; i++)
 688         {
 689             gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
 690         }
 691     }
 692
 693     /* Check that in the end we got the number of atoms we asked for */
 694     if (nrecv_local != ireq->n)
 695     {
 696         if (debug)
 697         {
 698             fprintf(debug, "Requested %d, received %d (tot recv %d)\n",
 699                     ireq->n, nrecv_local, nat_tot_specat-at_start);
 700             if (gmx_debug_at)
 701             {
 702                 for (i = 0; i < ireq->n; i++)
 703                 {
 704                     ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]);
 705                     fprintf(debug, " %s%d",
 706                             (ind >= 0) ? "" : "!",
 707                             ireq->ind[i]+1);
 708                 }
 709                 fprintf(debug, "\n");
 710             }
 711         }
 712         fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:",
 713                 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
 714         for (i = 0; i < ireq->n; i++)
 715         {
 716             if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0)
 717             {
 718                 fprintf(stderr, " %d", ireq->ind[i]+1);
 719             }
 720         }
 721         fprintf(stderr, "\n");
 722         gmx_fatal(FARGS, "DD cell %d %d %d could only obtain %d of the %d atoms that are connected via %ss from the neighboring cells. This probably means your %s lengths are too long compared to the domain decomposition cell size. Decrease the number of domain decomposition grid cells%s%s.",
 723                   dd->ci[XX], dd->ci[YY], dd->ci[ZZ],
 724                   nrecv_local, ireq->n, specat_type,
 725                   specat_type, add_err,
 726                   dd->bGridJump ? " or use the -rcon option of mdrun" : "");
 727     }
 728
 729     spac->at_start = at_start;
 730     spac->at_end   = nat_tot_specat;
 731
 732     if (debug)
 733     {
 734         fprintf(debug, "Done setup_specat_communication\n");
 735     }
 736
 737     return nat_tot_specat;
 738 }
 739
 740 static void walk_out(int con, int con_offset, int a, int offset, int nrec,
 741                      int ncon1, const t_iatom *ia1, const t_iatom *ia2,
 742                      const t_blocka *at2con,
 743                      const gmx_ga2la_t ga2la, gmx_bool bHomeConnect,
 744                      gmx_domdec_constraints_t *dc,
 745                      gmx_domdec_specat_comm_t *dcc,
 746                      t_ilist *il_local,
 747                      ind_req_t *ireq)
 748 {
 749     int            a1_gl, a2_gl, a_loc, i, coni, b;
 750     const t_iatom *iap;
 751
 752     if (dc->gc_req[con_offset+con] == 0)
 753     {
 754         /* Add this non-home constraint to the list */
 755         if (dc->ncon+1 > dc->con_nalloc)
 756         {
 757             dc->con_nalloc = over_alloc_large(dc->ncon+1);
 758             srenew(dc->con_gl, dc->con_nalloc);
 759             srenew(dc->con_nlocat, dc->con_nalloc);
 760         }
 761         dc->con_gl[dc->ncon]       = con_offset + con;
 762         dc->con_nlocat[dc->ncon]   = (bHomeConnect ? 1 : 0);
 763         dc->gc_req[con_offset+con] = 1;
 764         if (il_local->nr + 3 > il_local->nalloc)
 765         {
 766             il_local->nalloc = over_alloc_dd(il_local->nr+3);
 767             srenew(il_local->iatoms, il_local->nalloc);
 768         }
 769         iap = constr_iatomptr(ncon1, ia1, ia2, con);
 770         il_local->iatoms[il_local->nr++] = iap[0];
 771         a1_gl = offset + iap[1];
 772         a2_gl = offset + iap[2];
 773         /* The following indexing code can probably be optizimed */
 774         if (ga2la_get_home(ga2la, a1_gl, &a_loc))
 775         {
 776             il_local->iatoms[il_local->nr++] = a_loc;
 777         }
 778         else
 779         {
 780             /* We set this index later */
 781             il_local->iatoms[il_local->nr++] = -a1_gl - 1;
 782         }
 783         if (ga2la_get_home(ga2la, a2_gl, &a_loc))
 784         {
 785             il_local->iatoms[il_local->nr++] = a_loc;
 786         }
 787         else
 788         {
 789             /* We set this index later */
 790             il_local->iatoms[il_local->nr++] = -a2_gl - 1;
 791         }
 792         dc->ncon++;
 793     }
 794     /* Check to not ask for the same atom more than once */
 795     if (gmx_hash_get_minone(dc->ga2la, offset+a) == -1)
 796     {
 797         assert(dcc);
 798         /* Add this non-home atom to the list */
 799         if (ireq->n+1 > ireq->nalloc)
 800         {
 801             ireq->nalloc = over_alloc_large(ireq->n+1);
 802             srenew(ireq->ind, ireq->nalloc);
 803         }
 804         ireq->ind[ireq->n++] = offset + a;
 805         /* Temporarily mark with -2, we get the index later */
 806         gmx_hash_set(dc->ga2la, offset+a, -2);
 807     }
 808
 809     if (nrec > 0)
 810     {
 811         for (i = at2con->index[a]; i < at2con->index[a+1]; i++)
 812         {
 813             coni = at2con->a[i];
 814             if (coni != con)
 815             {
 816                 /* Walk further */
 817                 iap = constr_iatomptr(ncon1, ia1, ia2, coni);
 818                 if (a == iap[1])
 819                 {
 820                     b = iap[2];
 821                 }
 822                 else
 823                 {
 824                     b = iap[1];
 825                 }
 826                 if (!ga2la_get_home(ga2la, offset+b, &a_loc))
 827                 {
 828                     walk_out(coni, con_offset, b, offset, nrec-1,
 829                              ncon1, ia1, ia2, at2con,
 830                              ga2la, FALSE, dc, dcc, il_local, ireq);
 831                 }
 832             }
 833         }
 834     }
 835 }
 836
 837 static void atoms_to_settles(gmx_domdec_t *dd,
 838                              const gmx_mtop_t *mtop,
 839                              const int *cginfo,
 840                              const int **at2settle_mt,
 841                              int cg_start, int cg_end,
 842                              t_ilist *ils_local,
 843                              ind_req_t *ireq)
 844 {
 845     gmx_ga2la_t           ga2la;
 846     gmx_mtop_atomlookup_t alook;
 847     int                   settle;
 848     int                   nral, sa;
 849     int                   cg, a, a_gl, a_glsa, a_gls[3], a_locs[3];
 850     int                   mb, molnr, a_mol, offset;
 851     const gmx_molblock_t *molb;
 852     const t_iatom        *ia1;
 853     gmx_bool              a_home[3];
 854     int                   nlocal;
 855     gmx_bool              bAssign;
 856
 857     ga2la  = dd->ga2la;
 858
 859     alook = gmx_mtop_atomlookup_settle_init(mtop);
 860
 861     nral = NRAL(F_SETTLE);
 862
 863     for (cg = cg_start; cg < cg_end; cg++)
 864     {
 865         if (GET_CGINFO_SETTLE(cginfo[cg]))
 866         {
 867             for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
 868             {
 869                 a_gl = dd->gatindex[a];
 870
 871                 gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);
 872                 molb = &mtop->molblock[mb];
 873
 874                 settle = at2settle_mt[molb->type][a_mol];
 875
 876                 if (settle >= 0)
 877                 {
 878                     offset = a_gl - a_mol;
 879
 880                     ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms;
 881
 882                     bAssign = FALSE;
 883                     nlocal  = 0;
 884                     for (sa = 0; sa < nral; sa++)
 885                     {
 886                         a_glsa     = offset + ia1[settle*(1+nral)+1+sa];
 887                         a_gls[sa]  = a_glsa;
 888                         a_home[sa] = ga2la_get_home(ga2la, a_glsa, &a_locs[sa]);
 889                         if (a_home[sa])
 890                         {
 891                             if (nlocal == 0 && a_gl == a_glsa)
 892                             {
 893                                 bAssign = TRUE;
 894                             }
 895                             nlocal++;
 896                         }
 897                     }
 898
 899                     if (bAssign)
 900                     {
 901                         if (ils_local->nr+1+nral > ils_local->nalloc)
 902                         {
 903                             ils_local->nalloc = over_alloc_dd(ils_local->nr+1+nral);
 904                             srenew(ils_local->iatoms, ils_local->nalloc);
 905                         }
 906
 907                         ils_local->iatoms[ils_local->nr++] = ia1[settle*4];
 908
 909                         for (sa = 0; sa < nral; sa++)
 910                         {
 911                             if (ga2la_get_home(ga2la, a_gls[sa], &a_locs[sa]))
 912                             {
 913                                 ils_local->iatoms[ils_local->nr++] = a_locs[sa];
 914                             }
 915                             else
 916                             {
 917                                 ils_local->iatoms[ils_local->nr++] = -a_gls[sa] - 1;
 918                                 /* Add this non-home atom to the list */
 919                                 if (ireq->n+1 > ireq->nalloc)
 920                                 {
 921                                     ireq->nalloc = over_alloc_large(ireq->n+1);
 922                                     srenew(ireq->ind, ireq->nalloc);
 923                                 }
 924                                 ireq->ind[ireq->n++] = a_gls[sa];
 925                                 /* A check on double atom requests is
 926                                  * not required for settle.
 927                                  */
 928                             }
 929                         }
 930                     }
 931                 }
 932             }
 933         }
 934     }
 935
 936     gmx_mtop_atomlookup_destroy(alook);
 937 }
 938
 939 static void atoms_to_constraints(gmx_domdec_t *dd,
 940                                  const gmx_mtop_t *mtop,
 941                                  const int *cginfo,
 942                                  const t_blocka *at2con_mt, int nrec,
 943                                  t_ilist *ilc_local,
 944                                  ind_req_t *ireq)
 945 {
 946     const t_blocka           *at2con;
 947     gmx_ga2la_t               ga2la;
 948     gmx_mtop_atomlookup_t     alook;
 949     int                       ncon1;
 950     gmx_molblock_t           *molb;
 951     t_iatom                  *ia1, *ia2, *iap;
 952     int                       nhome, cg, a, a_gl, a_mol, a_loc, b_lo, offset, mb, molnr, b_mol, i, con, con_offset;
 953     gmx_domdec_constraints_t *dc;
 954     gmx_domdec_specat_comm_t *dcc;
 955
 956     dc  = dd->constraints;
 957     dcc = dd->constraint_comm;
 958
 959     ga2la  = dd->ga2la;
 960
 961     alook = gmx_mtop_atomlookup_init(mtop);
 962
 963     nhome = 0;
 964     for (cg = 0; cg < dd->ncg_home; cg++)
 965     {
 966         if (GET_CGINFO_CONSTR(cginfo[cg]))
 967         {
 968             for (a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
 969             {
 970                 a_gl = dd->gatindex[a];
 971
 972                 gmx_mtop_atomnr_to_molblock_ind(alook, a_gl, &mb, &molnr, &a_mol);
 973                 molb = &mtop->molblock[mb];
 974
 975                 ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/NRAL(F_SETTLE);
 976
 977                 ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms;
 978                 ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms;
 979
 980                 /* Calculate the global constraint number offset for the molecule.
 981                  * This is only required for the global index to make sure
 982                  * that we use each constraint only once.
 983                  */
 984                 con_offset =
 985                     dc->molb_con_offset[mb] + molnr*dc->molb_ncon_mol[mb];
 986
 987                 /* The global atom number offset for this molecule */
 988                 offset = a_gl - a_mol;
 989                 at2con = &at2con_mt[molb->type];
 990                 for (i = at2con->index[a_mol]; i < at2con->index[a_mol+1]; i++)
 991                 {
 992                     con = at2con->a[i];
 993                     iap = constr_iatomptr(ncon1, ia1, ia2, con);
 994                     if (a_mol == iap[1])
 995                     {
 996                         b_mol = iap[2];
 997                     }
 998                     else
 999                     {
1000                         b_mol = iap[1];
1001                     }
1002                     if (ga2la_get_home(ga2la, offset+b_mol, &a_loc))
1003                     {
1004                         /* Add this fully home constraint at the first atom */
1005                         if (a_mol < b_mol)
1006                         {
1007                             if (dc->ncon+1 > dc->con_nalloc)
1008                             {
1009                                 dc->con_nalloc = over_alloc_large(dc->ncon+1);
1010                                 srenew(dc->con_gl, dc->con_nalloc);
1011                                 srenew(dc->con_nlocat, dc->con_nalloc);
1012                             }
1013                             dc->con_gl[dc->ncon]     = con_offset + con;
1014                             dc->con_nlocat[dc->ncon] = 2;
1015                             if (ilc_local->nr + 3 > ilc_local->nalloc)
1016                             {
1017                                 ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3);
1018                                 srenew(ilc_local->iatoms, ilc_local->nalloc);
1019                             }
1020                             b_lo = a_loc;
1021                             ilc_local->iatoms[ilc_local->nr++] = iap[0];
1022                             ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? a    : b_lo);
1023                             ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? b_lo : a   );
1024                             dc->ncon++;
1025                             nhome++;
1026                         }
1027                     }
1028                     else
1029                     {
1030                         /* We need the nrec constraints coupled to this constraint,
1031                          * so we need to walk out of the home cell by nrec+1 atoms,
1032                          * since already atom bg is not locally present.
1033                          * Therefore we call walk_out with nrec recursions to go
1034                          * after this first call.
1035                          */
1036                         walk_out(con, con_offset, b_mol, offset, nrec,
1037                                  ncon1, ia1, ia2, at2con,
1038                                  dd->ga2la, TRUE, dc, dcc, ilc_local, ireq);
1039                     }
1040                 }
1041             }
1042         }
1043     }
1044
1045     gmx_mtop_atomlookup_destroy(alook);
1046
1047     if (debug)
1048     {
1049         fprintf(debug,
1050                 "Constraints: home %3d border %3d atoms: %3d\n",
1051                 nhome, dc->ncon-nhome,
1052                 dd->constraint_comm ? ireq->n : 0);
1053     }
1054 }
1055
1056 int dd_make_local_constraints(gmx_domdec_t *dd, int at_start,
1057                               const gmx_mtop_t *mtop,
1058                               const int *cginfo,
1059                               gmx_constr_t constr, int nrec,
1060                               t_ilist *il_local)
1061 {
1062     gmx_domdec_constraints_t *dc;
1063     t_ilist                  *ilc_local, *ils_local;
1064     ind_req_t                *ireq;
1065     const t_blocka           *at2con_mt;
1066     const int               **at2settle_mt;
1067     gmx_hash_t                ga2la_specat;
1068     int at_end, i, j;
1069     t_iatom                  *iap;
1070
1071     dc = dd->constraints;
1072
1073     ilc_local = &il_local[F_CONSTR];
1074     ils_local = &il_local[F_SETTLE];
1075
1076     dc->ncon      = 0;
1077     ilc_local->nr = 0;
1078     if (dd->constraint_comm)
1079     {
1080         at2con_mt = atom2constraints_moltype(constr);
1081         ireq      = &dd->constraint_comm->ireq[0];
1082         ireq->n   = 0;
1083     }
1084     else
1085     {
1086         at2con_mt = NULL;
1087         ireq      = NULL;
1088     }
1089
1090     if (dd->bInterCGsettles)
1091     {
1092         at2settle_mt  = atom2settle_moltype(constr);
1093         ils_local->nr = 0;
1094     }
1095     else
1096     {
1097         /* Settle works inside charge groups, we assigned them already */
1098         at2settle_mt = NULL;
1099     }
1100
1101     if (at2settle_mt == NULL)
1102     {
1103         atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1104                              ilc_local, ireq);
1105     }
1106     else
1107     {
1108         int t0_set;
1109         int thread;
1110
1111         /* Do the constraints, if present, on the first thread.
1112          * Do the settles on all other threads.
1113          */
1114         t0_set = ((at2con_mt != NULL && dc->nthread > 1) ? 1 : 0);
1115
1116 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
1117         for (thread = 0; thread < dc->nthread; thread++)
1118         {
1119             if (at2con_mt && thread == 0)
1120             {
1121                 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1122                                      ilc_local, ireq);
1123             }
1124
1125             if (thread >= t0_set)
1126             {
1127                 int        cg0, cg1;
1128                 t_ilist   *ilst;
1129                 ind_req_t *ireqt;
1130
1131                 /* Distribute the settle check+assignments over
1132                  * dc->nthread or dc->nthread-1 threads.
1133                  */
1134                 cg0 = (dd->ncg_home*(thread-t0_set  ))/(dc->nthread-t0_set);
1135                 cg1 = (dd->ncg_home*(thread-t0_set+1))/(dc->nthread-t0_set);
1136
1137                 if (thread == t0_set)
1138                 {
1139                     ilst = ils_local;
1140                 }
1141                 else
1142                 {
1143                     ilst = &dc->ils[thread];
1144                 }
1145                 ilst->nr = 0;
1146
1147                 ireqt = &dd->constraint_comm->ireq[thread];
1148                 if (thread > 0)
1149                 {
1150                     ireqt->n = 0;
1151                 }
1152
1153                 atoms_to_settles(dd, mtop, cginfo, at2settle_mt,
1154                                  cg0, cg1,
1155                                  ilst, ireqt);
1156             }
1157         }
1158
1159         /* Combine the generate settles and requested indices */
1160         for (thread = 1; thread < dc->nthread; thread++)
1161         {
1162             t_ilist   *ilst;
1163             ind_req_t *ireqt;
1164             int        ia;
1165
1166             if (thread > t0_set)
1167             {
1168                 ilst = &dc->ils[thread];
1169                 if (ils_local->nr + ilst->nr > ils_local->nalloc)
1170                 {
1171                     ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr);
1172                     srenew(ils_local->iatoms, ils_local->nalloc);
1173                 }
1174                 for (ia = 0; ia < ilst->nr; ia++)
1175                 {
1176                     ils_local->iatoms[ils_local->nr+ia] = ilst->iatoms[ia];
1177                 }
1178                 ils_local->nr += ilst->nr;
1179             }
1180
1181             ireqt = &dd->constraint_comm->ireq[thread];
1182             if (ireq->n+ireqt->n > ireq->nalloc)
1183             {
1184                 ireq->nalloc = over_alloc_large(ireq->n+ireqt->n);
1185                 srenew(ireq->ind, ireq->nalloc);
1186             }
1187             for (ia = 0; ia < ireqt->n; ia++)
1188             {
1189                 ireq->ind[ireq->n+ia] = ireqt->ind[ia];
1190             }
1191             ireq->n += ireqt->n;
1192         }
1193
1194         if (debug)
1195         {
1196             fprintf(debug, "Settles: total %3d\n", ils_local->nr/4);
1197         }
1198     }
1199
1200     if (dd->constraint_comm)
1201     {
1202         int nral1;
1203
1204         at_end =
1205             setup_specat_communication(dd, ireq, dd->constraint_comm,
1206                                        dd->constraints->ga2la,
1207                                        at_start, 2,
1208                                        "constraint", " or lincs-order");
1209
1210         /* Fill in the missing indices */
1211         ga2la_specat = dd->constraints->ga2la;
1212
1213         nral1 = 1 + NRAL(F_CONSTR);
1214         for (i = 0; i < ilc_local->nr; i += nral1)
1215         {
1216             iap = ilc_local->iatoms + i;
1217             for (j = 1; j < nral1; j++)
1218             {
1219                 if (iap[j] < 0)
1220                 {
1221                     iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1);
1222                 }
1223             }
1224         }
1225
1226         nral1 = 1 + NRAL(F_SETTLE);
1227         for (i = 0; i < ils_local->nr; i += nral1)
1228         {
1229             iap = ils_local->iatoms + i;
1230             for (j = 1; j < nral1; j++)
1231             {
1232                 if (iap[j] < 0)
1233                 {
1234                     iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1);
1235                 }
1236             }
1237         }
1238     }
1239     else
1240     {
1241         at_end = at_start;
1242     }
1243
1244     return at_end;
1245 }
1246
1247 int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
1248 {
1249     gmx_domdec_specat_comm_t *spac;
1250     ind_req_t                *ireq;
1251     gmx_hash_t                ga2la_specat;
1252     int  ftype, nral, i, j, gat, a;
1253     t_ilist                  *lilf;
1254     t_iatom                  *iatoms;
1255     int  at_end;
1256
1257     spac         = dd->vsite_comm;
1258     ireq         = &spac->ireq[0];
1259     ga2la_specat = dd->ga2la_vsite;
1260
1261     ireq->n = 0;
1262     /* Loop over all the home vsites */
1263     for (ftype = 0; ftype < F_NRE; ftype++)
1264     {
1265         if (interaction_function[ftype].flags & IF_VSITE)
1266         {
1267             nral = NRAL(ftype);
1268             lilf = &lil[ftype];
1269             for (i = 0; i < lilf->nr; i += 1+nral)
1270             {
1271                 iatoms = lilf->iatoms + i;
1272                 /* Check if we have the other atoms */
1273                 for (j = 1; j < 1+nral; j++)
1274                 {
1275                     if (iatoms[j] < 0)
1276                     {
1277                         /* This is not a home atom,
1278                          * we need to ask our neighbors.
1279                          */
1280                         a = -iatoms[j] - 1;
1281                         /* Check to not ask for the same atom more than once */
1282                         if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1)
1283                         {
1284                             /* Add this non-home atom to the list */
1285                             if (ireq->n+1 > ireq->nalloc)
1286                             {
1287                                 ireq->nalloc = over_alloc_large(ireq->n+1);
1288                                 srenew(ireq->ind, ireq->nalloc);
1289                             }
1290                             ireq->ind[ireq->n++] = a;
1291                             /* Temporarily mark with -2,
1292                              * we get the index later.
1293                              */
1294                             gmx_hash_set(ga2la_specat, a, -2);
1295                         }
1296                     }
1297                 }
1298             }
1299         }
1300     }
1301
1302     at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat,
1303                                         at_start, 1, "vsite", "");
1304
1305     /* Fill in the missing indices */
1306     for (ftype = 0; ftype < F_NRE; ftype++)
1307     {
1308         if (interaction_function[ftype].flags & IF_VSITE)
1309         {
1310             nral = NRAL(ftype);
1311             lilf = &lil[ftype];
1312             for (i = 0; i < lilf->nr; i += 1+nral)
1313             {
1314                 iatoms = lilf->iatoms + i;
1315                 for (j = 1; j < 1+nral; j++)
1316                 {
1317                     if (iatoms[j] < 0)
1318                     {
1319                         iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1);
1320                     }
1321                 }
1322             }
1323         }
1324     }
1325
1326     return at_end;
1327 }
1328
1329 static gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
1330 {
1331     gmx_domdec_specat_comm_t *spac;
1332
1333     snew(spac, 1);
1334     spac->nthread = nthread;
1335     snew(spac->ireq, spac->nthread);
1336
1337     return spac;
1338 }
1339
1340 void init_domdec_constraints(gmx_domdec_t *dd,
1341                              gmx_mtop_t   *mtop)
1342 {
1343     gmx_domdec_constraints_t *dc;
1344     gmx_molblock_t           *molb;
1345     int mb, ncon, c, a;
1346
1347     if (debug)
1348     {
1349         fprintf(debug, "Begin init_domdec_constraints\n");
1350     }
1351
1352     snew(dd->constraints, 1);
1353     dc = dd->constraints;
1354
1355     snew(dc->molb_con_offset, mtop->nmolblock);
1356     snew(dc->molb_ncon_mol, mtop->nmolblock);
1357
1358     ncon = 0;
1359     for (mb = 0; mb < mtop->nmolblock; mb++)
1360     {
1361         molb                    = &mtop->molblock[mb];
1362         dc->molb_con_offset[mb] = ncon;
1363         dc->molb_ncon_mol[mb]   =
1364             mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 +
1365             mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3;
1366         ncon += molb->nmol*dc->molb_ncon_mol[mb];
1367     }
1368
1369     if (ncon > 0)
1370     {
1371         snew(dc->gc_req, ncon);
1372         for (c = 0; c < ncon; c++)
1373         {
1374             dc->gc_req[c] = 0;
1375         }
1376     }
1377
1378     /* Use a hash table for the global to local index.
1379      * The number of keys is a rough estimate, it will be optimized later.
1380      */
1381     dc->ga2la = gmx_hash_init(min(mtop->natoms/20,
1382                                   mtop->natoms/(2*dd->nnodes)));
1383
1384     dc->nthread = gmx_omp_nthreads_get(emntDomdec);
1385     snew(dc->ils, dc->nthread);
1386
1387     dd->constraint_comm = specat_comm_init(dc->nthread);
1388 }
1389
1390 void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
1391 {
1392     int i;
1393     gmx_domdec_constraints_t *dc;
1394
1395     if (debug)
1396     {
1397         fprintf(debug, "Begin init_domdec_vsites\n");
1398     }
1399
1400     /* Use a hash table for the global to local index.
1401      * The number of keys is a rough estimate, it will be optimized later.
1402      */
1403     dd->ga2la_vsite = gmx_hash_init(min(n_intercg_vsite/20,
1404                                         n_intercg_vsite/(2*dd->nnodes)));
1405
1406     dd->vsite_comm = specat_comm_init(1);
1407 }