src/gromacs/mdlib/domdec_con.cpp

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