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 <assert.h>
  39
  40 #include <algorithm>
  41
  42 #include "gromacs/legacyheaders/constr.h"
  43 #include "gromacs/legacyheaders/domdec.h"
  44 #include "gromacs/legacyheaders/domdec_network.h"
  45 #include "gromacs/legacyheaders/gmx_ga2la.h"
  46 #include "gromacs/legacyheaders/gmx_hash.h"
  47 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
  48 #include "gromacs/legacyheaders/macros.h"
  49 #include "gromacs/legacyheaders/types/commrec.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/pbcutil/ishift.h"
  52 #include "gromacs/topology/mtop_util.h"
  53 #include "gromacs/utility/fatalerror.h"
  54 #include "gromacs/utility/gmxassert.h"
  55 #include "gromacs/utility/smalloc.h"
  56
  57 typedef struct {
  58     int  nsend;
  59     int *a;
  60     int  a_nalloc;
  61     int  nrecv;
  62 } gmx_specatsend_t;
  63
  64 typedef struct {
  65     int *ind;
  66     int  nalloc;
  67     int  n;
  68 } ind_req_t;
  69
  70 typedef struct gmx_domdec_specat_comm {
  71     /* The number of indices to receive during the setup */
  72     int              nreq[DIM][2][2];
  73     /* The atoms to send */
  74     gmx_specatsend_t spas[DIM][2];
  75     gmx_bool        *bSendAtom;
  76     int              bSendAtom_nalloc;
  77     /* Send buffers */
  78     int             *ibuf;
  79     int              ibuf_nalloc;
  80     rvec            *vbuf;
  81     int              vbuf_nalloc;
  82     rvec            *vbuf2;
  83     int              vbuf2_nalloc;
  84     /* The range in the local buffer(s) for received atoms */
  85     int              at_start;
  86     int              at_end;
  87
  88     /* The atom indices we need from the surrounding cells.
  89      * We can gather the indices over nthread threads.
  90      */
  91     int        nthread;
  92     ind_req_t *ireq;
  93 } gmx_domdec_specat_comm_t;
  94
  95 typedef struct gmx_domdec_constraints {
  96     int       *molb_con_offset;
  97     int       *molb_ncon_mol;
  98     /* The fully local and connected constraints */
  99     int        ncon;
 100     /* The global constraint number, only required for clearing gc_req */
 101     int       *con_gl;
 102     int       *con_nlocat;
 103     int        con_nalloc;
 104     /* Boolean that tells if a global constraint index has been requested */
 105     char      *gc_req;
 106     /* Global to local communicated constraint atom only index */
 107     gmx_hash_t ga2la;
 108
 109     /* Multi-threading stuff */
 110     int      nthread;
 111     t_ilist *ils;
 112 } gmx_domdec_constraints_t;
 113
 114
 115 static void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
 116                              rvec *f, rvec *fshift)
 117 {
 118     gmx_specatsend_t *spas;
 119     rvec             *vbuf;
 120     int               n, n0, n1, d, dim, dir, i;
 121     ivec              vis;
 122     int               is;
 123     gmx_bool          bPBC, bScrew;
 124
 125     n = spac->at_end;
 126     for (d = dd->ndim-1; d >= 0; d--)
 127     {
 128         dim = dd->dim[d];
 129         if (dd->nc[dim] > 2)
 130         {
 131             /* Pulse the grid forward and backward */
 132             spas = spac->spas[d];
 133             n0   = spas[0].nrecv;
 134             n1   = spas[1].nrecv;
 135             n   -= n1 + n0;
 136             vbuf = spac->vbuf;
 137             /* Send and receive the coordinates */
 138             dd_sendrecv2_rvec(dd, d,
 139                               f+n+n1, n0, vbuf, spas[0].nsend,
 140                               f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
 141             for (dir = 0; dir < 2; dir++)
 142             {
 143                 bPBC   = ((dir == 0 && dd->ci[dim] == 0) ||
 144                           (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
 145                 bScrew = (bPBC && dd->bScrewPBC && dim == XX);
 146
 147                 spas = &spac->spas[d][dir];
 148                 /* Sum the buffer into the required forces */
 149                 if (!bPBC || (!bScrew && fshift == NULL))
 150                 {
 151                     for (i = 0; i < spas->nsend; i++)
 152                     {
 153                         rvec_inc(f[spas->a[i]], *vbuf);
 154                         vbuf++;
 155                     }
 156                 }
 157                 else
 158                 {
 159                     clear_ivec(vis);
 160                     vis[dim] = (dir == 0 ? 1 : -1);
 161                     is       = IVEC2IS(vis);
 162                     if (!bScrew)
 163                     {
 164                         /* Sum and add to shift forces */
 165                         for (i = 0; i < spas->nsend; i++)
 166                         {
 167                             rvec_inc(f[spas->a[i]], *vbuf);
 168                             rvec_inc(fshift[is], *vbuf);
 169                             vbuf++;
 170                         }
 171                     }
 172                     else
 173                     {
 174                         /* Rotate the forces */
 175                         for (i = 0; i < spas->nsend; i++)
 176                         {
 177                             f[spas->a[i]][XX] += (*vbuf)[XX];
 178                             f[spas->a[i]][YY] -= (*vbuf)[YY];
 179                             f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
 180                             if (fshift)
 181                             {
 182                                 rvec_inc(fshift[is], *vbuf);
 183                             }
 184                             vbuf++;
 185                         }
 186                     }
 187                 }
 188             }
 189         }
 190         else
 191         {
 192             /* Two cells, so we only need to communicate one way */
 193             spas = &spac->spas[d][0];
 194             n   -= spas->nrecv;
 195             /* Send and receive the coordinates */
 196             dd_sendrecv_rvec(dd, d, dddirForward,
 197                              f+n, spas->nrecv, spac->vbuf, spas->nsend);
 198             /* Sum the buffer into the required forces */
 199             if (dd->bScrewPBC && dim == XX &&
 200                 (dd->ci[dim] == 0 ||
 201                  dd->ci[dim] == dd->nc[dim]-1))
 202             {
 203                 for (i = 0; i < spas->nsend; i++)
 204                 {
 205                     /* Rotate the force */
 206                     f[spas->a[i]][XX] += spac->vbuf[i][XX];
 207                     f[spas->a[i]][YY] -= spac->vbuf[i][YY];
 208                     f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
 209                 }
 210             }
 211             else
 212             {
 213                 for (i = 0; i < spas->nsend; i++)
 214                 {
 215                     rvec_inc(f[spas->a[i]], spac->vbuf[i]);
 216                 }
 217             }
 218         }
 219     }
 220 }
 221
 222 void dd_move_f_vsites(gmx_domdec_t *dd, rvec *f, rvec *fshift)
 223 {
 224     if (dd->vsite_comm)
 225     {
 226         dd_move_f_specat(dd, dd->vsite_comm, f, fshift);
 227     }
 228 }
 229
 230 void dd_clear_f_vsites(gmx_domdec_t *dd, rvec *f)
 231 {
 232     int i;
 233
 234     if (dd->vsite_comm)
 235     {
 236         for (i = dd->vsite_comm->at_start; i < dd->vsite_comm->at_end; i++)
 237         {
 238             clear_rvec(f[i]);
 239         }
 240     }
 241 }
 242
 243 static void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
 244                              matrix box,
 245                              rvec *x0,
 246                              rvec *x1, gmx_bool bX1IsCoord)
 247 {
 248     gmx_specatsend_t *spas;
 249     rvec             *x, *vbuf, *rbuf;
 250     int               nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
 251     gmx_bool          bPBC, bScrew = FALSE;
 252     rvec              shift = {0, 0, 0};
 253
 254     nvec = 1;
 255     if (x1 != NULL)
 256     {
 257         nvec++;
 258     }
 259
 260     n = spac->at_start;
 261     for (d = 0; d < dd->ndim; d++)
 262     {
 263         dim = dd->dim[d];
 264         if (dd->nc[dim] > 2)
 265         {
 266             /* Pulse the grid forward and backward */
 267             vbuf = spac->vbuf;
 268             for (dir = 0; dir < 2; dir++)
 269             {
 270                 if (dir == 0 && dd->ci[dim] == 0)
 271                 {
 272                     bPBC   = TRUE;
 273                     bScrew = (dd->bScrewPBC && dim == XX);
 274                     copy_rvec(box[dim], shift);
 275                 }
 276                 else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
 277                 {
 278                     bPBC   = TRUE;
 279                     bScrew = (dd->bScrewPBC && dim == XX);
 280                     for (i = 0; i < DIM; i++)
 281                     {
 282                         shift[i] = -box[dim][i];
 283                     }
 284                 }
 285                 else
 286                 {
 287                     bPBC   = FALSE;
 288                     bScrew = FALSE;
 289                 }
 290                 spas = &spac->spas[d][dir];
 291                 for (v = 0; v < nvec; v++)
 292                 {
 293                     x = (v == 0 ? x0 : x1);
 294                     /* Copy the required coordinates to the send buffer */
 295                     if (!bPBC || (v == 1 && !bX1IsCoord))
 296                     {
 297                         /* Only copy */
 298                         for (i = 0; i < spas->nsend; i++)
 299                         {
 300                             copy_rvec(x[spas->a[i]], *vbuf);
 301                             vbuf++;
 302                         }
 303                     }
 304                     else if (!bScrew)
 305                     {
 306                         /* Shift coordinates */
 307                         for (i = 0; i < spas->nsend; i++)
 308                         {
 309                             rvec_add(x[spas->a[i]], shift, *vbuf);
 310                             vbuf++;
 311                         }
 312                     }
 313                     else
 314                     {
 315                         /* Shift and rotate coordinates */
 316                         for (i = 0; i < spas->nsend; i++)
 317                         {
 318                             (*vbuf)[XX] =               x[spas->a[i]][XX] + shift[XX];
 319                             (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
 320                             (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
 321                             vbuf++;
 322                         }
 323                     }
 324                 }
 325             }
 326             /* Send and receive the coordinates */
 327             spas = spac->spas[d];
 328             ns0  = spas[0].nsend;
 329             nr0  = spas[0].nrecv;
 330             ns1  = spas[1].nsend;
 331             nr1  = spas[1].nrecv;
 332             if (nvec == 1)
 333             {
 334                 dd_sendrecv2_rvec(dd, d,
 335                                   spac->vbuf+ns0, ns1, x0+n, nr1,
 336                                   spac->vbuf, ns0, x0+n+nr1, nr0);
 337             }
 338             else
 339             {
 340                 /* Communicate both vectors in one buffer */
 341                 rbuf = spac->vbuf2;
 342                 dd_sendrecv2_rvec(dd, d,
 343                                   spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
 344                                   spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
 345                 /* Split the buffer into the two vectors */
 346                 nn = n;
 347                 for (dir = 1; dir >= 0; dir--)
 348                 {
 349                     nr = spas[dir].nrecv;
 350                     for (v = 0; v < 2; v++)
 351                     {
 352                         x = (v == 0 ? x0 : x1);
 353                         for (i = 0; i < nr; i++)
 354                         {
 355                             copy_rvec(*rbuf, x[nn+i]);
 356                             rbuf++;
 357                         }
 358                     }
 359                     nn += nr;
 360                 }
 361             }
 362             n += nr0 + nr1;
 363         }
 364         else
 365         {
 366             spas = &spac->spas[d][0];
 367             /* Copy the required coordinates to the send buffer */
 368             vbuf = spac->vbuf;
 369             for (v = 0; v < nvec; v++)
 370             {
 371                 x = (v == 0 ? x0 : x1);
 372                 if (dd->bScrewPBC && dim == XX &&
 373                     (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
 374                 {
 375                     /* Here we only perform the rotation, the rest of the pbc
 376                      * is handled in the constraint or viste routines.
 377                      */
 378                     for (i = 0; i < spas->nsend; i++)
 379                     {
 380                         (*vbuf)[XX] =               x[spas->a[i]][XX];
 381                         (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
 382                         (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
 383                         vbuf++;
 384                     }
 385                 }
 386                 else
 387                 {
 388                     for (i = 0; i < spas->nsend; i++)
 389                     {
 390                         copy_rvec(x[spas->a[i]], *vbuf);
 391                         vbuf++;
 392                     }
 393                 }
 394             }
 395             /* Send and receive the coordinates */
 396             if (nvec == 1)
 397             {
 398                 dd_sendrecv_rvec(dd, d, dddirBackward,
 399                                  spac->vbuf, spas->nsend, x0+n, spas->nrecv);
 400             }
 401             else
 402             {
 403                 /* Communicate both vectors in one buffer */
 404                 rbuf = spac->vbuf2;
 405                 dd_sendrecv_rvec(dd, d, dddirBackward,
 406                                  spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
 407                 /* Split the buffer into the two vectors */
 408                 nr = spas[0].nrecv;
 409                 for (v = 0; v < 2; v++)
 410                 {
 411                     x = (v == 0 ? x0 : x1);
 412                     for (i = 0; i < nr; i++)
 413                     {
 414                         copy_rvec(*rbuf, x[n+i]);
 415                         rbuf++;
 416                     }
 417                 }
 418             }
 419             n += spas->nrecv;
 420         }
 421     }
 422 }
 423
 424 void dd_move_x_constraints(gmx_domdec_t *dd, matrix box,
 425                            rvec *x0, rvec *x1, gmx_bool bX1IsCoord)
 426 {
 427     if (dd->constraint_comm)
 428     {
 429         dd_move_x_specat(dd, dd->constraint_comm, box, x0, x1, bX1IsCoord);
 430     }
 431 }
 432
 433 void dd_move_x_vsites(gmx_domdec_t *dd, matrix box, rvec *x)
 434 {
 435     if (dd->vsite_comm)
 436     {
 437         dd_move_x_specat(dd, dd->vsite_comm, box, x, NULL, FALSE);
 438     }
 439 }
 440
 441 int *dd_constraints_nlocalatoms(gmx_domdec_t *dd)
 442 {
 443     if (dd->constraints)
 444     {
 445         return dd->constraints->con_nlocat;
 446     }
 447     else
 448     {
 449         return NULL;
 450     }
 451 }
 452
 453 void dd_clear_local_constraint_indices(gmx_domdec_t *dd)
 454 {
 455     gmx_domdec_constraints_t *dc;
 456     int i;
 457
 458     dc = dd->constraints;
 459
 460     for (i = 0; i < dc->ncon; i++)
 461     {
 462         dc->gc_req[dc->con_gl[i]] = 0;
 463     }
 464
 465     if (dd->constraint_comm)
 466     {
 467         gmx_hash_clear_and_optimize(dc->ga2la);
 468     }
 469 }
 470
 471 void dd_clear_local_vsite_indices(gmx_domdec_t *dd)
 472 {
 473     if (dd->vsite_comm)
 474     {
 475         gmx_hash_clear_and_optimize(dd->ga2la_vsite);
 476     }
 477 }
 478
 479 static int setup_specat_communication(gmx_domdec_t             *dd,
 480                                       ind_req_t                *ireq,
 481                                       gmx_domdec_specat_comm_t *spac,
 482                                       gmx_hash_t                ga2la_specat,
 483                                       int                       at_start,
 484                                       int                       vbuf_fac,
 485                                       const char               *specat_type,
 486                                       const char               *add_err)
 487 {
 488     int               nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
 489     int               d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
 490     int               nat_tot_specat, nat_tot_prev, nalloc_old;
 491     gmx_bool          bPBC;
 492     gmx_specatsend_t *spas;
 493
 494     if (debug)
 495     {
 496         fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
 497     }
 498
 499     /* nsend[0]: the number of atoms requested by this node only,
 500      *           we communicate this for more efficients checks
 501      * nsend[1]: the total number of requested atoms
 502      */
 503     nsend[0] = ireq->n;
 504     nsend[1] = nsend[0];
 505     nlast    = nsend[1];
 506     for (d = dd->ndim-1; d >= 0; d--)
 507     {
 508         /* Pulse the grid forward and backward */
 509         dim  = dd->dim[d];
 510         bPBC = (dim < dd->npbcdim);
 511         if (dd->nc[dim] == 2)
 512         {
 513             /* Only 2 cells, so we only need to communicate once */
 514             ndir = 1;
 515         }
 516         else
 517         {
 518             ndir = 2;
 519         }
 520         for (dir = 0; dir < ndir; dir++)
 521         {
 522             if (!bPBC &&
 523                 dd->nc[dim] > 2 &&
 524                 ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
 525                  (dir == 1 && dd->ci[dim] == 0)))
 526             {
 527                 /* No pbc: the fist/last cell should not request atoms */
 528                 nsend_ptr = nsend_zero;
 529             }
 530             else
 531             {
 532                 nsend_ptr = nsend;
 533             }
 534             /* Communicate the number of indices */
 535             dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
 536                             nsend_ptr, 2, spac->nreq[d][dir], 2);
 537             nr = spac->nreq[d][dir][1];
 538             if (nlast+nr > ireq->nalloc)
 539             {
 540                 ireq->nalloc = over_alloc_dd(nlast+nr);
 541                 srenew(ireq->ind, ireq->nalloc);
 542             }
 543             /* Communicate the indices */
 544             dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
 545                             ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
 546             nlast += nr;
 547         }
 548         nsend[1] = nlast;
 549     }
 550     if (debug)
 551     {
 552         fprintf(debug, "Communicated the counts\n");
 553     }
 554
 555     /* Search for the requested atoms and communicate the indices we have */
 556     nat_tot_specat = at_start;
 557     nrecv_local    = 0;
 558     for (d = 0; d < dd->ndim; d++)
 559     {
 560         /* 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     // This code should not be called unless this condition is true,
1077     // because that's the only time init_domdec_constraints is
1078     // called...
1079     GMX_RELEASE_ASSERT(dd->bInterCGcons || dd->bInterCGsettles, "dd_make_local_constraints called when there are no local constraints");
1080     // ... and init_domdec_constraints always sets
1081     // dd->constraint_comm...
1082     GMX_RELEASE_ASSERT(dd->constraint_comm, "Invalid use of dd_make_local_constraints before construction of constraint_comm");
1083     // ... which static analysis needs to be reassured about, because
1084     // otherwise, when dd->bInterCGsettles is
1085     // true. dd->constraint_comm is unilaterally dereferenced before
1086     // the call to atoms_to_settles.
1087
1088     dc = dd->constraints;
1089
1090     ilc_local = &il_local[F_CONSTR];
1091     ils_local = &il_local[F_SETTLE];
1092
1093     dc->ncon      = 0;
1094     ilc_local->nr = 0;
1095     if (dd->constraint_comm)
1096     {
1097         at2con_mt = atom2constraints_moltype(constr);
1098         ireq      = &dd->constraint_comm->ireq[0];
1099         ireq->n   = 0;
1100     }
1101     else
1102     {
1103         // Currently unreachable
1104         at2con_mt = NULL;
1105         ireq      = NULL;
1106     }
1107
1108     if (dd->bInterCGsettles)
1109     {
1110         at2settle_mt  = atom2settle_moltype(constr);
1111         ils_local->nr = 0;
1112     }
1113     else
1114     {
1115         /* Settle works inside charge groups, we assigned them already */
1116         at2settle_mt = NULL;
1117     }
1118
1119     if (at2settle_mt == NULL)
1120     {
1121         atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1122                              ilc_local, ireq);
1123     }
1124     else
1125     {
1126         int t0_set;
1127         int thread;
1128
1129         /* Do the constraints, if present, on the first thread.
1130          * Do the settles on all other threads.
1131          */
1132         t0_set = ((at2con_mt != NULL && dc->nthread > 1) ? 1 : 0);
1133
1134 #pragma omp parallel for num_threads(dc->nthread) schedule(static)
1135         for (thread = 0; thread < dc->nthread; thread++)
1136         {
1137             if (at2con_mt && thread == 0)
1138             {
1139                 atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
1140                                      ilc_local, ireq);
1141             }
1142
1143             if (thread >= t0_set)
1144             {
1145                 int        cg0, cg1;
1146                 t_ilist   *ilst;
1147                 ind_req_t *ireqt;
1148
1149                 /* Distribute the settle check+assignments over
1150                  * dc->nthread or dc->nthread-1 threads.
1151                  */
1152                 cg0 = (dd->ncg_home*(thread-t0_set  ))/(dc->nthread-t0_set);
1153                 cg1 = (dd->ncg_home*(thread-t0_set+1))/(dc->nthread-t0_set);
1154
1155                 if (thread == t0_set)
1156                 {
1157                     ilst = ils_local;
1158                 }
1159                 else
1160                 {
1161                     ilst = &dc->ils[thread];
1162                 }
1163                 ilst->nr = 0;
1164
1165                 ireqt = &dd->constraint_comm->ireq[thread];
1166                 if (thread > 0)
1167                 {
1168                     ireqt->n = 0;
1169                 }
1170
1171                 atoms_to_settles(dd, mtop, cginfo, at2settle_mt,
1172                                  cg0, cg1,
1173                                  ilst, ireqt);
1174             }
1175         }
1176
1177         /* Combine the generate settles and requested indices */
1178         for (thread = 1; thread < dc->nthread; thread++)
1179         {
1180             t_ilist   *ilst;
1181             ind_req_t *ireqt;
1182             int        ia;
1183
1184             if (thread > t0_set)
1185             {
1186                 ilst = &dc->ils[thread];
1187                 if (ils_local->nr + ilst->nr > ils_local->nalloc)
1188                 {
1189                     ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr);
1190                     srenew(ils_local->iatoms, ils_local->nalloc);
1191                 }
1192                 for (ia = 0; ia < ilst->nr; ia++)
1193                 {
1194                     ils_local->iatoms[ils_local->nr+ia] = ilst->iatoms[ia];
1195                 }
1196                 ils_local->nr += ilst->nr;
1197             }
1198
1199             ireqt = &dd->constraint_comm->ireq[thread];
1200             if (ireq->n+ireqt->n > ireq->nalloc)
1201             {
1202                 ireq->nalloc = over_alloc_large(ireq->n+ireqt->n);
1203                 srenew(ireq->ind, ireq->nalloc);
1204             }
1205             for (ia = 0; ia < ireqt->n; ia++)
1206             {
1207                 ireq->ind[ireq->n+ia] = ireqt->ind[ia];
1208             }
1209             ireq->n += ireqt->n;
1210         }
1211
1212         if (debug)
1213         {
1214             fprintf(debug, "Settles: total %3d\n", ils_local->nr/4);
1215         }
1216     }
1217
1218     if (dd->constraint_comm)
1219     {
1220         int nral1;
1221
1222         at_end =
1223             setup_specat_communication(dd, ireq, dd->constraint_comm,
1224                                        dd->constraints->ga2la,
1225                                        at_start, 2,
1226                                        "constraint", " or lincs-order");
1227
1228         /* Fill in the missing indices */
1229         ga2la_specat = dd->constraints->ga2la;
1230
1231         nral1 = 1 + NRAL(F_CONSTR);
1232         for (i = 0; i < ilc_local->nr; i += nral1)
1233         {
1234             iap = ilc_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         nral1 = 1 + NRAL(F_SETTLE);
1245         for (i = 0; i < ils_local->nr; i += nral1)
1246         {
1247             iap = ils_local->iatoms + i;
1248             for (j = 1; j < nral1; j++)
1249             {
1250                 if (iap[j] < 0)
1251                 {
1252                     iap[j] = gmx_hash_get_minone(ga2la_specat, -iap[j]-1);
1253                 }
1254             }
1255         }
1256     }
1257     else
1258     {
1259         // Currently unreachable
1260         at_end = at_start;
1261     }
1262
1263     return at_end;
1264 }
1265
1266 int dd_make_local_vsites(gmx_domdec_t *dd, int at_start, t_ilist *lil)
1267 {
1268     gmx_domdec_specat_comm_t *spac;
1269     ind_req_t                *ireq;
1270     gmx_hash_t                ga2la_specat;
1271     int  ftype, nral, i, j, a;
1272     t_ilist                  *lilf;
1273     t_iatom                  *iatoms;
1274     int  at_end;
1275
1276     spac         = dd->vsite_comm;
1277     ireq         = &spac->ireq[0];
1278     ga2la_specat = dd->ga2la_vsite;
1279
1280     ireq->n = 0;
1281     /* Loop over all the home vsites */
1282     for (ftype = 0; ftype < F_NRE; ftype++)
1283     {
1284         if (interaction_function[ftype].flags & IF_VSITE)
1285         {
1286             nral = NRAL(ftype);
1287             lilf = &lil[ftype];
1288             for (i = 0; i < lilf->nr; i += 1+nral)
1289             {
1290                 iatoms = lilf->iatoms + i;
1291                 /* Check if we have the other atoms */
1292                 for (j = 1; j < 1+nral; j++)
1293                 {
1294                     if (iatoms[j] < 0)
1295                     {
1296                         /* This is not a home atom,
1297                          * we need to ask our neighbors.
1298                          */
1299                         a = -iatoms[j] - 1;
1300                         /* Check to not ask for the same atom more than once */
1301                         if (gmx_hash_get_minone(dd->ga2la_vsite, a) == -1)
1302                         {
1303                             /* Add this non-home atom to the list */
1304                             if (ireq->n+1 > ireq->nalloc)
1305                             {
1306                                 ireq->nalloc = over_alloc_large(ireq->n+1);
1307                                 srenew(ireq->ind, ireq->nalloc);
1308                             }
1309                             ireq->ind[ireq->n++] = a;
1310                             /* Temporarily mark with -2,
1311                              * we get the index later.
1312                              */
1313                             gmx_hash_set(ga2la_specat, a, -2);
1314                         }
1315                     }
1316                 }
1317             }
1318         }
1319     }
1320
1321     at_end = setup_specat_communication(dd, ireq, dd->vsite_comm, ga2la_specat,
1322                                         at_start, 1, "vsite", "");
1323
1324     /* Fill in the missing indices */
1325     for (ftype = 0; ftype < F_NRE; ftype++)
1326     {
1327         if (interaction_function[ftype].flags & IF_VSITE)
1328         {
1329             nral = NRAL(ftype);
1330             lilf = &lil[ftype];
1331             for (i = 0; i < lilf->nr; i += 1+nral)
1332             {
1333                 iatoms = lilf->iatoms + i;
1334                 for (j = 1; j < 1+nral; j++)
1335                 {
1336                     if (iatoms[j] < 0)
1337                     {
1338                         iatoms[j] = gmx_hash_get_minone(ga2la_specat, -iatoms[j]-1);
1339                     }
1340                 }
1341             }
1342         }
1343     }
1344
1345     return at_end;
1346 }
1347
1348 static gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
1349 {
1350     gmx_domdec_specat_comm_t *spac;
1351
1352     snew(spac, 1);
1353     spac->nthread = nthread;
1354     snew(spac->ireq, spac->nthread);
1355
1356     return spac;
1357 }
1358
1359 void init_domdec_constraints(gmx_domdec_t *dd,
1360                              gmx_mtop_t   *mtop)
1361 {
1362     gmx_domdec_constraints_t *dc;
1363     gmx_molblock_t           *molb;
1364     int mb, ncon, c;
1365
1366     if (debug)
1367     {
1368         fprintf(debug, "Begin init_domdec_constraints\n");
1369     }
1370
1371     snew(dd->constraints, 1);
1372     dc = dd->constraints;
1373
1374     snew(dc->molb_con_offset, mtop->nmolblock);
1375     snew(dc->molb_ncon_mol, mtop->nmolblock);
1376
1377     ncon = 0;
1378     for (mb = 0; mb < mtop->nmolblock; mb++)
1379     {
1380         molb                    = &mtop->molblock[mb];
1381         dc->molb_con_offset[mb] = ncon;
1382         dc->molb_ncon_mol[mb]   =
1383             mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 +
1384             mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3;
1385         ncon += molb->nmol*dc->molb_ncon_mol[mb];
1386     }
1387
1388     if (ncon > 0)
1389     {
1390         snew(dc->gc_req, ncon);
1391         for (c = 0; c < ncon; c++)
1392         {
1393             dc->gc_req[c] = 0;
1394         }
1395     }
1396
1397     /* Use a hash table for the global to local index.
1398      * The number of keys is a rough estimate, it will be optimized later.
1399      */
1400     dc->ga2la = gmx_hash_init(std::min(mtop->natoms/20,
1401                                        mtop->natoms/(2*dd->nnodes)));
1402
1403     dc->nthread = gmx_omp_nthreads_get(emntDomdec);
1404     snew(dc->ils, dc->nthread);
1405
1406     dd->constraint_comm = specat_comm_init(dc->nthread);
1407 }
1408
1409 void init_domdec_vsites(gmx_domdec_t *dd, int n_intercg_vsite)
1410 {
1411     if (debug)
1412     {
1413         fprintf(debug, "Begin init_domdec_vsites\n");
1414     }
1415
1416     /* Use a hash table for the global to local index.
1417      * The number of keys is a rough estimate, it will be optimized later.
1418      */
1419     dd->ga2la_vsite = gmx_hash_init(std::min(n_intercg_vsite/20,
1420                                              n_intercg_vsite/(2*dd->nnodes)));
1421
1422     dd->vsite_comm = specat_comm_init(1);
1423 }