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 by the GROMACS development team.
   5  * Copyright (c) 2012,2013,2014,2015,2017 by the GROMACS development team.
   6  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 /*! \internal \file
  39  *
  40  * \brief This file implements functions for domdec to use
  41  * while managing inter-atomic constraints.
  42  *
  43  * \author Berk Hess <hess@kth.se>
  44  * \ingroup module_domdec
  45  */
  46
  47 #include "gmxpre.h"
  48
  49 #include "domdec_specatomcomm.h"
  50
  51 #include <cassert>
  52
  53 #include <algorithm>
  54
  55 #include "gromacs/domdec/dlb.h"
  56 #include "gromacs/domdec/domdec.h"
  57 #include "gromacs/domdec/domdec_network.h"
  58 #include "gromacs/domdec/domdec_struct.h"
  59 #include "gromacs/domdec/ga2la.h"
  60 #include "gromacs/domdec/hashedmap.h"
  61 #include "gromacs/domdec/partition.h"
  62 #include "gromacs/math/vec.h"
  63 #include "gromacs/mdtypes/commrec.h"
  64 #include "gromacs/pbcutil/ishift.h"
  65 #include "gromacs/utility/fatalerror.h"
  66 #include "gromacs/utility/gmxassert.h"
  67
  68 void dd_move_f_specat(const gmx_domdec_t* dd, gmx_domdec_specat_comm_t* spac, rvec* f, rvec* fshift)
  69 {
  70     gmx_specatsend_t* spas;
  71     rvec*             vbuf;
  72     int               n, n0, n1, dim, dir;
  73     ivec              vis;
  74     int               is;
  75     gmx_bool          bPBC, bScrew;
  76
  77     n = spac->at_end;
  78     for (int d = dd->ndim - 1; d >= 0; d--)
  79     {
  80         dim = dd->dim[d];
  81         if (dd->numCells[dim] > 2)
  82         {
  83             /* Pulse the grid forward and backward */
  84             spas = spac->spas[d];
  85             n0   = spas[0].nrecv;
  86             n1   = spas[1].nrecv;
  87             n -= n1 + n0;
  88             vbuf = as_rvec_array(spac->vbuf.data());
  89             /* Send and receive the coordinates */
  90             dd_sendrecv2_rvec(dd, d, f + n + n1, n0, vbuf, spas[0].a.size(), f + n, n1,
  91                               vbuf + spas[0].a.size(), spas[1].a.size());
  92             for (dir = 0; dir < 2; dir++)
  93             {
  94                 bPBC   = ((dir == 0 && dd->ci[dim] == 0)
  95                         || (dir == 1 && dd->ci[dim] == dd->numCells[dim] - 1));
  96                 bScrew = (bPBC && dd->unitCellInfo.haveScrewPBC && dim == XX);
  97
  98                 spas = &spac->spas[d][dir];
  99                 /* Sum the buffer into the required forces */
 100                 if (!bPBC || (!bScrew && fshift == nullptr))
 101                 {
 102                     for (int a : spas->a)
 103                     {
 104                         rvec_inc(f[a], *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 (int a : spas->a)
 117                         {
 118                             rvec_inc(f[a], *vbuf);
 119                             rvec_inc(fshift[is], *vbuf);
 120                             vbuf++;
 121                         }
 122                     }
 123                     else
 124                     {
 125                         /* Rotate the forces */
 126                         for (int a : spas->a)
 127                         {
 128                             f[a][XX] += (*vbuf)[XX];
 129                             f[a][YY] -= (*vbuf)[YY];
 130                             f[a][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             ddSendrecv(dd, d, dddirForward, f + n, spas->nrecv, as_rvec_array(spac->vbuf.data()),
 148                        spas->a.size());
 149             /* Sum the buffer into the required forces */
 150             if (dd->unitCellInfo.haveScrewPBC && dim == XX
 151                 && (dd->ci[dim] == 0 || dd->ci[dim] == dd->numCells[dim] - 1))
 152             {
 153                 int i = 0;
 154                 for (int a : spas->a)
 155                 {
 156                     /* Rotate the force */
 157                     f[a][XX] += spac->vbuf[i][XX];
 158                     f[a][YY] -= spac->vbuf[i][YY];
 159                     f[a][ZZ] -= spac->vbuf[i][ZZ];
 160                     i++;
 161                 }
 162             }
 163             else
 164             {
 165                 int i = 0;
 166                 for (int a : spas->a)
 167                 {
 168                     rvec_inc(f[a], spac->vbuf[i]);
 169                     i++;
 170                 }
 171             }
 172         }
 173     }
 174 }
 175
 176 void dd_move_x_specat(const gmx_domdec_t*       dd,
 177                       gmx_domdec_specat_comm_t* spac,
 178                       const matrix              box,
 179                       rvec*                     x0,
 180                       rvec*                     x1,
 181                       gmx_bool                  bX1IsCoord)
 182 {
 183     gmx_specatsend_t* spas;
 184     int               nvec, v, n, nn, ns0, ns1, nr0, nr1, nr, d, dim, dir, i;
 185     gmx_bool          bPBC, bScrew = FALSE;
 186     rvec              shift = { 0, 0, 0 };
 187
 188     nvec = 1;
 189     if (x1 != nullptr)
 190     {
 191         nvec++;
 192     }
 193
 194     n = spac->at_start;
 195     for (d = 0; d < dd->ndim; d++)
 196     {
 197         dim = dd->dim[d];
 198         if (dd->numCells[dim] > 2)
 199         {
 200             /* Pulse the grid forward and backward */
 201             rvec* vbuf = as_rvec_array(spac->vbuf.data());
 202             for (dir = 0; dir < 2; dir++)
 203             {
 204                 if (dir == 0 && dd->ci[dim] == 0)
 205                 {
 206                     bPBC   = TRUE;
 207                     bScrew = (dd->unitCellInfo.haveScrewPBC && dim == XX);
 208                     copy_rvec(box[dim], shift);
 209                 }
 210                 else if (dir == 1 && dd->ci[dim] == dd->numCells[dim] - 1)
 211                 {
 212                     bPBC   = TRUE;
 213                     bScrew = (dd->unitCellInfo.haveScrewPBC && dim == XX);
 214                     for (i = 0; i < DIM; i++)
 215                     {
 216                         shift[i] = -box[dim][i];
 217                     }
 218                 }
 219                 else
 220                 {
 221                     bPBC   = FALSE;
 222                     bScrew = FALSE;
 223                 }
 224                 spas = &spac->spas[d][dir];
 225                 for (v = 0; v < nvec; v++)
 226                 {
 227                     rvec* x = (v == 0 ? x0 : x1);
 228                     /* Copy the required coordinates to the send buffer */
 229                     if (!bPBC || (v == 1 && !bX1IsCoord))
 230                     {
 231                         /* Only copy */
 232                         for (int a : spas->a)
 233                         {
 234                             copy_rvec(x[a], *vbuf);
 235                             vbuf++;
 236                         }
 237                     }
 238                     else if (!bScrew)
 239                     {
 240                         /* Shift coordinates */
 241                         for (int a : spas->a)
 242                         {
 243                             rvec_add(x[a], shift, *vbuf);
 244                             vbuf++;
 245                         }
 246                     }
 247                     else
 248                     {
 249                         /* Shift and rotate coordinates */
 250                         for (int a : spas->a)
 251                         {
 252                             (*vbuf)[XX] = x[a][XX] + shift[XX];
 253                             (*vbuf)[YY] = box[YY][YY] - x[a][YY] + shift[YY];
 254                             (*vbuf)[ZZ] = box[ZZ][ZZ] - x[a][ZZ] + shift[ZZ];
 255                             vbuf++;
 256                         }
 257                     }
 258                 }
 259             }
 260             /* Send and receive the coordinates */
 261             spas = spac->spas[d];
 262             ns0  = spas[0].a.size();
 263             nr0  = spas[0].nrecv;
 264             ns1  = spas[1].a.size();
 265             nr1  = spas[1].nrecv;
 266             if (nvec == 1)
 267             {
 268                 rvec* vbuf = as_rvec_array(spac->vbuf.data());
 269                 dd_sendrecv2_rvec(dd, d, vbuf + ns0, ns1, x0 + n, nr1, vbuf, ns0, x0 + n + nr1, nr0);
 270             }
 271             else
 272             {
 273                 rvec* vbuf = as_rvec_array(spac->vbuf.data());
 274                 /* Communicate both vectors in one buffer */
 275                 rvec* rbuf = as_rvec_array(spac->vbuf2.data());
 276                 dd_sendrecv2_rvec(dd, d, vbuf + 2 * ns0, 2 * ns1, rbuf, 2 * nr1, vbuf, 2 * ns0,
 277                                   rbuf + 2 * nr1, 2 * nr0);
 278                 /* Split the buffer into the two vectors */
 279                 nn = n;
 280                 for (dir = 1; dir >= 0; dir--)
 281                 {
 282                     nr = spas[dir].nrecv;
 283                     for (v = 0; v < 2; v++)
 284                     {
 285                         rvec* x = (v == 0 ? x0 : x1);
 286                         for (i = 0; i < nr; i++)
 287                         {
 288                             copy_rvec(*rbuf, x[nn + i]);
 289                             rbuf++;
 290                         }
 291                     }
 292                     nn += nr;
 293                 }
 294             }
 295             n += nr0 + nr1;
 296         }
 297         else
 298         {
 299             spas = &spac->spas[d][0];
 300             /* Copy the required coordinates to the send buffer */
 301             rvec* vbuf = as_rvec_array(spac->vbuf.data());
 302             for (v = 0; v < nvec; v++)
 303             {
 304                 rvec* x = (v == 0 ? x0 : x1);
 305                 if (dd->unitCellInfo.haveScrewPBC && dim == XX
 306                     && (dd->ci[XX] == 0 || dd->ci[XX] == dd->numCells[XX] - 1))
 307                 {
 308                     /* Here we only perform the rotation, the rest of the pbc
 309                      * is handled in the constraint or viste routines.
 310                      */
 311                     for (int a : spas->a)
 312                     {
 313                         (*vbuf)[XX] = x[a][XX];
 314                         (*vbuf)[YY] = box[YY][YY] - x[a][YY];
 315                         (*vbuf)[ZZ] = box[ZZ][ZZ] - x[a][ZZ];
 316                         vbuf++;
 317                     }
 318                 }
 319                 else
 320                 {
 321                     for (int a : spas->a)
 322                     {
 323                         copy_rvec(x[a], *vbuf);
 324                         vbuf++;
 325                     }
 326                 }
 327             }
 328             /* Send and receive the coordinates */
 329             if (nvec == 1)
 330             {
 331                 rvec* vbuf = as_rvec_array(spac->vbuf.data());
 332                 ddSendrecv(dd, d, dddirBackward, vbuf, spas->a.size(), x0 + n, spas->nrecv);
 333             }
 334             else
 335             {
 336                 rvec* vbuf = as_rvec_array(spac->vbuf.data());
 337                 /* Communicate both vectors in one buffer */
 338                 rvec* rbuf = as_rvec_array(spac->vbuf2.data());
 339                 ddSendrecv(dd, d, dddirBackward, vbuf, 2 * spas->a.size(), rbuf, 2 * spas->nrecv);
 340                 /* Split the buffer into the two vectors */
 341                 nr = spas[0].nrecv;
 342                 for (v = 0; v < 2; v++)
 343                 {
 344                     rvec* x = (v == 0 ? x0 : x1);
 345                     for (i = 0; i < nr; i++)
 346                     {
 347                         copy_rvec(*rbuf, x[n + i]);
 348                         rbuf++;
 349                     }
 350                 }
 351             }
 352             n += spas->nrecv;
 353         }
 354     }
 355 }
 356
 357 int setup_specat_communication(gmx_domdec_t*             dd,
 358                                std::vector<int>*         ireq,
 359                                gmx_domdec_specat_comm_t* spac,
 360                                gmx::HashedMap<int>*      ga2la_specat,
 361                                int                       at_start,
 362                                int                       vbuf_fac,
 363                                const char*               specat_type,
 364                                const char*               add_err)
 365 {
 366     int               nsend[2], nlast, nsend_zero[2] = { 0, 0 }, *nsend_ptr;
 367     int               dim, ndir, nr, ns, nrecv_local, n0, start, buf[2];
 368     int               nat_tot_specat, nat_tot_prev;
 369     gmx_bool          bPBC;
 370     gmx_specatsend_t* spas;
 371
 372     if (debug)
 373     {
 374         fprintf(debug, "Begin setup_specat_communication for %s\n", specat_type);
 375     }
 376
 377     /* nsend[0]: the number of atoms requested by this node only,
 378      *           we communicate this for more efficients checks
 379      * nsend[1]: the total number of requested atoms
 380      */
 381     const int numRequested = ireq->size();
 382     nsend[0]               = ireq->size();
 383     nsend[1]               = nsend[0];
 384     nlast                  = nsend[1];
 385     for (int d = dd->ndim - 1; d >= 0; d--)
 386     {
 387         /* Pulse the grid forward and backward */
 388         dim  = dd->dim[d];
 389         bPBC = (dim < dd->unitCellInfo.npbcdim);
 390         if (dd->numCells[dim] == 2)
 391         {
 392             /* Only 2 cells, so we only need to communicate once */
 393             ndir = 1;
 394         }
 395         else
 396         {
 397             ndir = 2;
 398         }
 399         for (int dir = 0; dir < ndir; dir++)
 400         {
 401             if (!bPBC && dd->numCells[dim] > 2
 402                 && ((dir == 0 && dd->ci[dim] == dd->numCells[dim] - 1) || (dir == 1 && dd->ci[dim] == 0)))
 403             {
 404                 /* No pbc: the fist/last cell should not request atoms */
 405                 nsend_ptr = nsend_zero;
 406             }
 407             else
 408             {
 409                 nsend_ptr = nsend;
 410             }
 411             /* Communicate the number of indices */
 412             ddSendrecv(dd, d, dir == 0 ? dddirForward : dddirBackward, nsend_ptr, 2,
 413                        spac->nreq[d][dir], 2);
 414             nr = spac->nreq[d][dir][1];
 415             ireq->resize(nlast + nr);
 416             /* Communicate the indices */
 417             ddSendrecv(dd, d, dir == 0 ? dddirForward : dddirBackward, ireq->data(), nsend_ptr[1],
 418                        ireq->data() + nlast, nr);
 419             nlast += nr;
 420         }
 421         nsend[1] = nlast;
 422     }
 423     if (debug)
 424     {
 425         fprintf(debug, "Communicated the counts\n");
 426     }
 427
 428     /* Search for the requested atoms and communicate the indices we have */
 429     nat_tot_specat = at_start;
 430     nrecv_local    = 0;
 431     for (int d = 0; d < dd->ndim; d++)
 432     {
 433         /* Pulse the grid forward and backward */
 434         if (dd->dim[d] >= dd->unitCellInfo.npbcdim || dd->numCells[dd->dim[d]] > 2)
 435         {
 436             ndir = 2;
 437         }
 438         else
 439         {
 440             ndir = 1;
 441         }
 442         nat_tot_prev = nat_tot_specat;
 443         for (int dir = ndir - 1; dir >= 0; dir--)
 444         {
 445             /* To avoid cost of clearing by resize(), we only increase size */
 446             if (static_cast<size_t>(nat_tot_specat) > spac->sendAtom.size())
 447             {
 448                 /* Note: resize initializes new elements to false, which is actually needed here */
 449                 spac->sendAtom.resize(nat_tot_specat);
 450             }
 451             spas = &spac->spas[d][dir];
 452             n0   = spac->nreq[d][dir][0];
 453             nr   = spac->nreq[d][dir][1];
 454             if (debug)
 455             {
 456                 fprintf(debug, "dim=%d, dir=%d, searching for %d atoms\n", d, dir, nr);
 457             }
 458             start = nlast - nr;
 459             spas->a.clear();
 460             spac->ibuf.clear();
 461             nsend[0] = 0;
 462             for (int i = 0; i < nr; i++)
 463             {
 464                 const int indr = (*ireq)[start + i];
 465                 int       ind;
 466                 /* Check if this is a home atom and if so ind will be set */
 467                 if (const int* homeIndex = dd->ga2la->findHome(indr))
 468                 {
 469                     ind = *homeIndex;
 470                 }
 471                 else
 472                 {
 473                     /* Search in the communicated atoms */
 474                     if (const int* a = ga2la_specat->find(indr))
 475                     {
 476                         ind = *a;
 477                     }
 478                     else
 479                     {
 480                         ind = -1;
 481                     }
 482                 }
 483                 if (ind >= 0)
 484                 {
 485                     if (i < n0 || !spac->sendAtom[ind])
 486                     {
 487                         /* Store the local index so we know which coordinates
 488                          * to send out later.
 489                          */
 490                         spas->a.push_back(ind);
 491                         spac->sendAtom[ind] = true;
 492                         /* Store the global index so we can send it now */
 493                         spac->ibuf.push_back(indr);
 494                         if (i < n0)
 495                         {
 496                             nsend[0]++;
 497                         }
 498                     }
 499                 }
 500             }
 501             nlast = start;
 502             /* Clear the local flags */
 503             for (int a : spas->a)
 504             {
 505                 spac->sendAtom[a] = false;
 506             }
 507             /* Send and receive the number of indices to communicate */
 508             nsend[1] = spas->a.size();
 509             ddSendrecv(dd, d, dir == 0 ? dddirBackward : dddirForward, nsend, 2, buf, 2);
 510             if (debug)
 511             {
 512                 fprintf(debug,
 513                         "Send to rank %d, %d (%d) indices, "
 514                         "receive from rank %d, %d (%d) indices\n",
 515                         dd->neighbor[d][1 - dir], nsend[1], nsend[0], dd->neighbor[d][dir], buf[1],
 516                         buf[0]);
 517                 if (gmx_debug_at)
 518                 {
 519                     for (int i : spac->ibuf)
 520                     {
 521                         fprintf(debug, " %d", i + 1);
 522                     }
 523                     fprintf(debug, "\n");
 524                 }
 525             }
 526             nrecv_local += buf[0];
 527             spas->nrecv = buf[1];
 528             dd->globalAtomIndices.resize(nat_tot_specat + spas->nrecv);
 529             /* Send and receive the indices */
 530             ddSendrecv(dd, d, dir == 0 ? dddirBackward : dddirForward, spac->ibuf.data(),
 531                        spac->ibuf.size(), dd->globalAtomIndices.data() + nat_tot_specat, spas->nrecv);
 532             nat_tot_specat += spas->nrecv;
 533         }
 534
 535         /* Increase the x/f communication buffer sizes, when necessary */
 536         ns = spac->spas[d][0].a.size();
 537         nr = spac->spas[d][0].nrecv;
 538         if (ndir == 2)
 539         {
 540             ns += spac->spas[d][1].a.size();
 541             nr += spac->spas[d][1].nrecv;
 542         }
 543         if (vbuf_fac * ns > gmx::index(spac->vbuf.size()))
 544         {
 545             spac->vbuf.resize(vbuf_fac * ns);
 546         }
 547         if (vbuf_fac == 2 && vbuf_fac * nr > gmx::index(spac->vbuf2.size()))
 548         {
 549             spac->vbuf2.resize(vbuf_fac * nr);
 550         }
 551
 552         /* Make a global to local index for the communication atoms */
 553         for (int i = nat_tot_prev; i < nat_tot_specat; i++)
 554         {
 555             ga2la_specat->insert_or_assign(dd->globalAtomIndices[i], i);
 556         }
 557     }
 558
 559     /* Check that in the end we got the number of atoms we asked for */
 560     if (nrecv_local != numRequested)
 561     {
 562         if (debug)
 563         {
 564             fprintf(debug, "Requested %d, received %d (tot recv %d)\n", numRequested, nrecv_local,
 565                     nat_tot_specat - at_start);
 566             if (gmx_debug_at)
 567             {
 568                 for (int i = 0; i < numRequested; i++)
 569                 {
 570                     const int* ind = ga2la_specat->find((*ireq)[i]);
 571                     fprintf(debug, " %s%d", ind ? "" : "!", (*ireq)[i] + 1);
 572                 }
 573                 fprintf(debug, "\n");
 574             }
 575         }
 576         fprintf(stderr, "\nDD cell %d %d %d: Neighboring cells do not have atoms:", dd->ci[XX],
 577                 dd->ci[YY], dd->ci[ZZ]);
 578         for (int i = 0; i < numRequested; i++)
 579         {
 580             if (!ga2la_specat->find((*ireq)[i]))
 581             {
 582                 fprintf(stderr, " %d", (*ireq)[i] + 1);
 583             }
 584         }
 585         fprintf(stderr, "\n");
 586         gmx_fatal(FARGS,
 587                   "DD cell %d %d %d could only obtain %d of the %d atoms that are connected via "
 588                   "%ss from the neighboring cells. This probably means your %s lengths are too "
 589                   "long compared to the domain decomposition cell size. Decrease the number of "
 590                   "domain decomposition grid cells%s%s.",
 591                   dd->ci[XX], dd->ci[YY], dd->ci[ZZ], nrecv_local, numRequested, specat_type,
 592                   specat_type, add_err, dd_dlb_is_on(dd) ? " or use the -rcon option of mdrun" : "");
 593     }
 594
 595     spac->at_start = at_start;
 596     spac->at_end   = nat_tot_specat;
 597
 598     if (debug)
 599     {
 600         fprintf(debug, "Done setup_specat_communication\n");
 601     }
 602
 603     return nat_tot_specat;
 604 }