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