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