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