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