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