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