src/gromacs/domdec/domdec_specatomcomm.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,2015, 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 /*! \internal \file
  37  *
  38  * \brief This file implements functions for domdec to use
  39  * while managing inter-atomic constraints.
  40  *
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_domdec
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "domdec_specatomcomm.h"
  48
  49 #include <assert.h>
  50
  51 #include <algorithm>
  52
  53 #include "gromacs/domdec/domdec.h"
  54 #include "gromacs/domdec/domdec_network.h"
  55 #include "gromacs/legacyheaders/gmx_ga2la.h"
  56 #include "gromacs/legacyheaders/gmx_hash.h"
  57 #include "gromacs/legacyheaders/gmx_omp_nthreads.h"
  58 #include "gromacs/legacyheaders/macros.h"
  59 #include "gromacs/legacyheaders/types/commrec.h"
  60 #include "gromacs/math/vec.h"
  61 #include "gromacs/pbcutil/ishift.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/gmxassert.h"
  64 #include "gromacs/utility/smalloc.h"
  65
  66 void dd_move_f_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
  67                       rvec *f, rvec *fshift)
  68 {
  69     gmx_specatsend_t *spas;
  70     rvec             *vbuf;
  71     int               n, n0, n1, d, dim, dir, i;
  72     ivec              vis;
  73     int               is;
  74     gmx_bool          bPBC, bScrew;
  75
  76     n = spac->at_end;
  77     for (d = dd->ndim-1; d >= 0; d--)
  78     {
  79         dim = dd->dim[d];
  80         if (dd->nc[dim] > 2)
  81         {
  82             /* Pulse the grid forward and backward */
  83             spas = spac->spas[d];
  84             n0   = spas[0].nrecv;
  85             n1   = spas[1].nrecv;
  86             n   -= n1 + n0;
  87             vbuf = spac->vbuf;
  88             /* Send and receive the coordinates */
  89             dd_sendrecv2_rvec(dd, d,
  90                               f+n+n1, n0, vbuf, spas[0].nsend,
  91                               f+n, n1, vbuf+spas[0].nsend, spas[1].nsend);
  92             for (dir = 0; dir < 2; dir++)
  93             {
  94                 bPBC   = ((dir == 0 && dd->ci[dim] == 0) ||
  95                           (dir == 1 && dd->ci[dim] == dd->nc[dim]-1));
  96                 bScrew = (bPBC && dd->bScrewPBC && dim == XX);
  97
  98                 spas = &spac->spas[d][dir];
  99                 /* Sum the buffer into the required forces */
 100                 if (!bPBC || (!bScrew && fshift == NULL))
 101                 {
 102                     for (i = 0; i < spas->nsend; i++)
 103                     {
 104                         rvec_inc(f[spas->a[i]], *vbuf);
 105                         vbuf++;
 106                     }
 107                 }
 108                 else
 109                 {
 110                     clear_ivec(vis);
 111                     vis[dim] = (dir == 0 ? 1 : -1);
 112                     is       = IVEC2IS(vis);
 113                     if (!bScrew)
 114                     {
 115                         /* Sum and add to shift forces */
 116                         for (i = 0; i < spas->nsend; i++)
 117                         {
 118                             rvec_inc(f[spas->a[i]], *vbuf);
 119                             rvec_inc(fshift[is], *vbuf);
 120                             vbuf++;
 121                         }
 122                     }
 123                     else
 124                     {
 125                         /* Rotate the forces */
 126                         for (i = 0; i < spas->nsend; i++)
 127                         {
 128                             f[spas->a[i]][XX] += (*vbuf)[XX];
 129                             f[spas->a[i]][YY] -= (*vbuf)[YY];
 130                             f[spas->a[i]][ZZ] -= (*vbuf)[ZZ];
 131                             if (fshift)
 132                             {
 133                                 rvec_inc(fshift[is], *vbuf);
 134                             }
 135                             vbuf++;
 136                         }
 137                     }
 138                 }
 139             }
 140         }
 141         else
 142         {
 143             /* Two cells, so we only need to communicate one way */
 144             spas = &spac->spas[d][0];
 145             n   -= spas->nrecv;
 146             /* Send and receive the coordinates */
 147             dd_sendrecv_rvec(dd, d, dddirForward,
 148                              f+n, spas->nrecv, spac->vbuf, spas->nsend);
 149             /* Sum the buffer into the required forces */
 150             if (dd->bScrewPBC && dim == XX &&
 151                 (dd->ci[dim] == 0 ||
 152                  dd->ci[dim] == dd->nc[dim]-1))
 153             {
 154                 for (i = 0; i < spas->nsend; i++)
 155                 {
 156                     /* Rotate the force */
 157                     f[spas->a[i]][XX] += spac->vbuf[i][XX];
 158                     f[spas->a[i]][YY] -= spac->vbuf[i][YY];
 159                     f[spas->a[i]][ZZ] -= spac->vbuf[i][ZZ];
 160                 }
 161             }
 162             else
 163             {
 164                 for (i = 0; i < spas->nsend; i++)
 165                 {
 166                     rvec_inc(f[spas->a[i]], spac->vbuf[i]);
 167                 }
 168             }
 169         }
 170     }
 171 }
 172
 173 void dd_move_x_specat(gmx_domdec_t *dd, gmx_domdec_specat_comm_t *spac,
 174                       matrix box,
 175                       rvec *x0,
 176                       rvec *x1, gmx_bool bX1IsCoord)
 177 {
 178     gmx_specatsend_t *spas;
 179     rvec             *x, *vbuf, *rbuf;
 180     int               nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
 181     gmx_bool          bPBC, bScrew = FALSE;
 182     rvec              shift = {0, 0, 0};
 183
 184     nvec = 1;
 185     if (x1 != NULL)
 186     {
 187         nvec++;
 188     }
 189
 190     n = spac->at_start;
 191     for (d = 0; d < dd->ndim; d++)
 192     {
 193         dim = dd->dim[d];
 194         if (dd->nc[dim] > 2)
 195         {
 196             /* Pulse the grid forward and backward */
 197             vbuf = spac->vbuf;
 198             for (dir = 0; dir < 2; dir++)
 199             {
 200                 if (dir == 0 && dd->ci[dim] == 0)
 201                 {
 202                     bPBC   = TRUE;
 203                     bScrew = (dd->bScrewPBC && dim == XX);
 204                     copy_rvec(box[dim], shift);
 205                 }
 206                 else if (dir == 1 && dd->ci[dim] == dd->nc[dim]-1)
 207                 {
 208                     bPBC   = TRUE;
 209                     bScrew = (dd->bScrewPBC && dim == XX);
 210                     for (i = 0; i < DIM; i++)
 211                     {
 212                         shift[i] = -box[dim][i];
 213                     }
 214                 }
 215                 else
 216                 {
 217                     bPBC   = FALSE;
 218                     bScrew = FALSE;
 219                 }
 220                 spas = &spac->spas[d][dir];
 221                 for (v = 0; v < nvec; v++)
 222                 {
 223                     x = (v == 0 ? x0 : x1);
 224                     /* Copy the required coordinates to the send buffer */
 225                     if (!bPBC || (v == 1 && !bX1IsCoord))
 226                     {
 227                         /* Only copy */
 228                         for (i = 0; i < spas->nsend; i++)
 229                         {
 230                             copy_rvec(x[spas->a[i]], *vbuf);
 231                             vbuf++;
 232                         }
 233                     }
 234                     else if (!bScrew)
 235                     {
 236                         /* Shift coordinates */
 237                         for (i = 0; i < spas->nsend; i++)
 238                         {
 239                             rvec_add(x[spas->a[i]], shift, *vbuf);
 240                             vbuf++;
 241                         }
 242                     }
 243                     else
 244                     {
 245                         /* Shift and rotate coordinates */
 246                         for (i = 0; i < spas->nsend; i++)
 247                         {
 248                             (*vbuf)[XX] =               x[spas->a[i]][XX] + shift[XX];
 249                             (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY] + shift[YY];
 250                             (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ] + shift[ZZ];
 251                             vbuf++;
 252                         }
 253                     }
 254                 }
 255             }
 256             /* Send and receive the coordinates */
 257             spas = spac->spas[d];
 258             ns0  = spas[0].nsend;
 259             nr0  = spas[0].nrecv;
 260             ns1  = spas[1].nsend;
 261             nr1  = spas[1].nrecv;
 262             if (nvec == 1)
 263             {
 264                 dd_sendrecv2_rvec(dd, d,
 265                                   spac->vbuf+ns0, ns1, x0+n, nr1,
 266                                   spac->vbuf, ns0, x0+n+nr1, nr0);
 267             }
 268             else
 269             {
 270                 /* Communicate both vectors in one buffer */
 271                 rbuf = spac->vbuf2;
 272                 dd_sendrecv2_rvec(dd, d,
 273                                   spac->vbuf+2*ns0, 2*ns1, rbuf, 2*nr1,
 274                                   spac->vbuf, 2*ns0, rbuf+2*nr1, 2*nr0);
 275                 /* Split the buffer into the two vectors */
 276                 nn = n;
 277                 for (dir = 1; dir >= 0; dir--)
 278                 {
 279                     nr = spas[dir].nrecv;
 280                     for (v = 0; v < 2; v++)
 281                     {
 282                         x = (v == 0 ? x0 : x1);
 283                         for (i = 0; i < nr; i++)
 284                         {
 285                             copy_rvec(*rbuf, x[nn+i]);
 286                             rbuf++;
 287                         }
 288                     }
 289                     nn += nr;
 290                 }
 291             }
 292             n += nr0 + nr1;
 293         }
 294         else
 295         {
 296             spas = &spac->spas[d][0];
 297             /* Copy the required coordinates to the send buffer */
 298             vbuf = spac->vbuf;
 299             for (v = 0; v < nvec; v++)
 300             {
 301                 x = (v == 0 ? x0 : x1);
 302                 if (dd->bScrewPBC && dim == XX &&
 303                     (dd->ci[XX] == 0 || dd->ci[XX] == dd->nc[XX]-1))
 304                 {
 305                     /* Here we only perform the rotation, the rest of the pbc
 306                      * is handled in the constraint or viste routines.
 307                      */
 308                     for (i = 0; i < spas->nsend; i++)
 309                     {
 310                         (*vbuf)[XX] =               x[spas->a[i]][XX];
 311                         (*vbuf)[YY] = box[YY][YY] - x[spas->a[i]][YY];
 312                         (*vbuf)[ZZ] = box[ZZ][ZZ] - x[spas->a[i]][ZZ];
 313                         vbuf++;
 314                     }
 315                 }
 316                 else
 317                 {
 318                     for (i = 0; i < spas->nsend; i++)
 319                     {
 320                         copy_rvec(x[spas->a[i]], *vbuf);
 321                         vbuf++;
 322                     }
 323                 }
 324             }
 325             /* Send and receive the coordinates */
 326             if (nvec == 1)
 327             {
 328                 dd_sendrecv_rvec(dd, d, dddirBackward,
 329                                  spac->vbuf, spas->nsend, x0+n, spas->nrecv);
 330             }
 331             else
 332             {
 333                 /* Communicate both vectors in one buffer */
 334                 rbuf = spac->vbuf2;
 335                 dd_sendrecv_rvec(dd, d, dddirBackward,
 336                                  spac->vbuf, 2*spas->nsend, rbuf, 2*spas->nrecv);
 337                 /* Split the buffer into the two vectors */
 338                 nr = spas[0].nrecv;
 339                 for (v = 0; v < 2; v++)
 340                 {
 341                     x = (v == 0 ? x0 : x1);
 342                     for (i = 0; i < nr; i++)
 343                     {
 344                         copy_rvec(*rbuf, x[n+i]);
 345                         rbuf++;
 346                     }
 347                 }
 348             }
 349             n += spas->nrecv;
 350         }
 351     }
 352 }
 353
 354 int setup_specat_communication(gmx_domdec_t             *dd,
 355                                ind_req_t                *ireq,
 356                                gmx_domdec_specat_comm_t *spac,
 357                                gmx_hash_t                ga2la_specat,
 358                                int                       at_start,
 359                                int                       vbuf_fac,
 360                                const char               *specat_type,
 361                                const char               *add_err)
 362 {
 363     int               nsend[2], nlast, nsend_zero[2] = {0, 0}, *nsend_ptr;
 364     int               d, dim, ndir, dir, nr, ns, i, nrecv_local, n0, start, indr, ind, buf[2];
 365     int               nat_tot_specat, nat_tot_prev, nalloc_old;
 366     gmx_bool          bPBC;
 367     gmx_specatsend_t *spas;
 368
 369     if (debug)
 370     {
 371         fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
 372     }
 373
 374     /* nsend[0]: the number of atoms requested by this node only,
 375      *           we communicate this for more efficients checks
 376      * nsend[1]: the total number of requested atoms
 377      */
 378     nsend[0] = ireq->n;
 379     nsend[1] = nsend[0];
 380     nlast    = nsend[1];
 381     for (d = dd->ndim-1; d >= 0; d--)
 382     {
 383         /* Pulse the grid forward and backward */
 384         dim  = dd->dim[d];
 385         bPBC = (dim < dd->npbcdim);
 386         if (dd->nc[dim] == 2)
 387         {
 388             /* Only 2 cells, so we only need to communicate once */
 389             ndir = 1;
 390         }
 391         else
 392         {
 393             ndir = 2;
 394         }
 395         for (dir = 0; dir < ndir; dir++)
 396         {
 397             if (!bPBC &&
 398                 dd->nc[dim] > 2 &&
 399                 ((dir == 0 && dd->ci[dim] == dd->nc[dim] - 1) ||
 400                  (dir == 1 && dd->ci[dim] == 0)))
 401             {
 402                 /* No pbc: the fist/last cell should not request atoms */
 403                 nsend_ptr = nsend_zero;
 404             }
 405             else
 406             {
 407                 nsend_ptr = nsend;
 408             }
 409             /* Communicate the number of indices */
 410             dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
 411                             nsend_ptr, 2, spac->nreq[d][dir], 2);
 412             nr = spac->nreq[d][dir][1];
 413             if (nlast+nr > ireq->nalloc)
 414             {
 415                 ireq->nalloc = over_alloc_dd(nlast+nr);
 416                 srenew(ireq->ind, ireq->nalloc);
 417             }
 418             /* Communicate the indices */
 419             dd_sendrecv_int(dd, d, dir == 0 ? dddirForward : dddirBackward,
 420                             ireq->ind, nsend_ptr[1], ireq->ind+nlast, nr);
 421             nlast += nr;
 422         }
 423         nsend[1] = nlast;
 424     }
 425     if (debug)
 426     {
 427         fprintf(debug, "Communicated the counts\n");
 428     }
 429
 430     /* Search for the requested atoms and communicate the indices we have */
 431     nat_tot_specat = at_start;
 432     nrecv_local    = 0;
 433     for (d = 0; d < dd->ndim; d++)
 434     {
 435         /* Pulse the grid forward and backward */
 436         if (dd->dim[d] >= dd->npbcdim || dd->nc[dd->dim[d]] > 2)
 437         {
 438             ndir = 2;
 439         }
 440         else
 441         {
 442             ndir = 1;
 443         }
 444         nat_tot_prev = nat_tot_specat;
 445         for (dir = ndir-1; dir >= 0; dir--)
 446         {
 447             if (nat_tot_specat > spac->bSendAtom_nalloc)
 448             {
 449                 nalloc_old             = spac->bSendAtom_nalloc;
 450                 spac->bSendAtom_nalloc = over_alloc_dd(nat_tot_specat);
 451                 srenew(spac->bSendAtom, spac->bSendAtom_nalloc);
 452                 for (i = nalloc_old; i < spac->bSendAtom_nalloc; i++)
 453                 {
 454                     spac->bSendAtom[i] = FALSE;
 455                 }
 456             }
 457             spas = &spac->spas[d][dir];
 458             n0   = spac->nreq[d][dir][0];
 459             nr   = spac->nreq[d][dir][1];
 460             if (debug)
 461             {
 462                 fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n",
 463                         d, dir, nr);
 464             }
 465             start       = nlast - nr;
 466             spas->nsend = 0;
 467             nsend[0]    = 0;
 468             for (i = 0; i < nr; i++)
 469             {
 470                 indr = ireq->ind[start+i];
 471                 ind  = -1;
 472                 /* Check if this is a home atom and if so ind will be set */
 473                 if (!ga2la_get_home(dd->ga2la, indr, &ind))
 474                 {
 475                     /* Search in the communicated atoms */
 476                     ind = gmx_hash_get_minone(ga2la_specat, indr);
 477                 }
 478                 if (ind >= 0)
 479                 {
 480                     if (i < n0 || !spac->bSendAtom[ind])
 481                     {
 482                         if (spas->nsend+1 > spas->a_nalloc)
 483                         {
 484                             spas->a_nalloc = over_alloc_large(spas->nsend+1);
 485                             srenew(spas->a, spas->a_nalloc);
 486                         }
 487                         /* Store the local index so we know which coordinates
 488                          * to send out later.
 489                          */
 490                         spas->a[spas->nsend] = ind;
 491                         spac->bSendAtom[ind] = TRUE;
 492                         if (spas->nsend+1 > spac->ibuf_nalloc)
 493                         {
 494                             spac->ibuf_nalloc = over_alloc_large(spas->nsend+1);
 495                             srenew(spac->ibuf, spac->ibuf_nalloc);
 496                         }
 497                         /* Store the global index so we can send it now */
 498                         spac->ibuf[spas->nsend] = indr;
 499                         if (i < n0)
 500                         {
 501                             nsend[0]++;
 502                         }
 503                         spas->nsend++;
 504                     }
 505                 }
 506             }
 507             nlast = start;
 508             /* Clear the local flags */
 509             for (i = 0; i < spas->nsend; i++)
 510             {
 511                 spac->bSendAtom[spas->a[i]] = FALSE;
 512             }
 513             /* Send and receive the number of indices to communicate */
 514             nsend[1] = spas->nsend;
 515             dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
 516                             nsend, 2, buf, 2);
 517             if (debug)
 518             {
 519                 fprintf(debug, "Send to rank %d, %d (%d) indices, "
 520                         "receive from rank %d, %d (%d) indices\n",
 521                         dd->neighbor[d][1-dir], nsend[1], nsend[0],
 522                         dd->neighbor[d][dir], buf[1], buf[0]);
 523                 if (gmx_debug_at)
 524                 {
 525                     for (i = 0; i < spas->nsend; i++)
 526                     {
 527                         fprintf(debug, " %d", spac->ibuf[i]+1);
 528                     }
 529                     fprintf(debug, "\n");
 530                 }
 531             }
 532             nrecv_local += buf[0];
 533             spas->nrecv  = buf[1];
 534             if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
 535             {
 536                 dd->gatindex_nalloc =
 537                     over_alloc_dd(nat_tot_specat + spas->nrecv);
 538                 srenew(dd->gatindex, dd->gatindex_nalloc);
 539             }
 540             /* Send and receive the indices */
 541             dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
 542                             spac->ibuf, spas->nsend,
 543                             dd->gatindex+nat_tot_specat, spas->nrecv);
 544             nat_tot_specat += spas->nrecv;
 545         }
 546
 547         /* Allocate the x/f communication buffers */
 548         ns = spac->spas[d][0].nsend;
 549         nr = spac->spas[d][0].nrecv;
 550         if (ndir == 2)
 551         {
 552             ns += spac->spas[d][1].nsend;
 553             nr += spac->spas[d][1].nrecv;
 554         }
 555         if (vbuf_fac*ns > spac->vbuf_nalloc)
 556         {
 557             spac->vbuf_nalloc = over_alloc_dd(vbuf_fac*ns);
 558             srenew(spac->vbuf, spac->vbuf_nalloc);
 559         }
 560         if (vbuf_fac == 2 && vbuf_fac*nr > spac->vbuf2_nalloc)
 561         {
 562             spac->vbuf2_nalloc = over_alloc_dd(vbuf_fac*nr);
 563             srenew(spac->vbuf2, spac->vbuf2_nalloc);
 564         }
 565
 566         /* Make a global to local index for the communication atoms */
 567         for (i = nat_tot_prev; i < nat_tot_specat; i++)
 568         {
 569             gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
 570         }
 571     }
 572
 573     /* Check that in the end we got the number of atoms we asked for */
 574     if (nrecv_local != ireq->n)
 575     {
 576         if (debug)
 577         {
 578             fprintf(debug, "Requested %d, received %d (tot recv %d)\n",
 579                     ireq->n, nrecv_local, nat_tot_specat-at_start);
 580             if (gmx_debug_at)
 581             {
 582                 for (i = 0; i < ireq->n; i++)
 583                 {
 584                     ind = gmx_hash_get_minone(ga2la_specat, ireq->ind[i]);
 585                     fprintf(debug, " %s%d",
 586                             (ind >= 0) ? "" : "!",
 587                             ireq->ind[i]+1);
 588                 }
 589                 fprintf(debug, "\n");
 590             }
 591         }
 592         fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:",
 593                 dd->ci[XX], dd->ci[YY], dd->ci[ZZ]);
 594         for (i = 0; i < ireq->n; i++)
 595         {
 596             if (gmx_hash_get_minone(ga2la_specat, ireq->ind[i]) < 0)
 597             {
 598                 fprintf(stderr, " %d", ireq->ind[i]+1);
 599             }
 600         }
 601         fprintf(stderr, "\n");
 602         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.",
 603                   dd->ci[XX], dd->ci[YY], dd->ci[ZZ],
 604                   nrecv_local, ireq->n, specat_type,
 605                   specat_type, add_err,
 606                   dd_dlb_is_on(dd) ? " or use the -rcon option of mdrun" : "");
 607     }
 608
 609     spac->at_start = at_start;
 610     spac->at_end   = nat_tot_specat;
 611
 612     if (debug)
 613     {
 614         fprintf(debug, "Done setup_specat_communication\n");
 615     }
 616
 617     return nat_tot_specat;
 618 }
 619
 620 gmx_domdec_specat_comm_t *specat_comm_init(int nthread)
 621 {
 622     gmx_domdec_specat_comm_t *spac;
 623
 624     snew(spac, 1);
 625     spac->nthread = nthread;
 626     snew(spac->ireq, spac->nthread);
 627
 628     return spac;
 629 }