src/mdlib/domdec_con.c

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