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