src/gromacs/domdec/domdec_topology.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2006 - 2014, The GROMACS development team.
   5  * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36
  37 /*! \internal \file
  38  *
  39  * \brief This file defines functions used by the domdec module
  40  * while managing the construction, use and error checking for
  41  * topologies local to a DD rank.
  42  *
  43  * \author Berk Hess <hess@kth.se>
  44  * \ingroup module_domdec
  45  */
  46
  47 #include "gmxpre.h"
  48
  49 #include <cassert>
  50 #include <cstdlib>
  51 #include <cstring>
  52
  53 #include <algorithm>
  54 #include <memory>
  55 #include <optional>
  56 #include <string>
  57
  58 #include "gromacs/domdec/domdec.h"
  59 #include "gromacs/domdec/domdec_network.h"
  60 #include "gromacs/domdec/ga2la.h"
  61 #include "gromacs/domdec/localtopologychecker.h"
  62 #include "gromacs/domdec/options.h"
  63 #include "gromacs/domdec/reversetopology.h"
  64 #include "gromacs/gmxlib/network.h"
  65 #include "gromacs/math/vec.h"
  66 #include "gromacs/mdlib/forcerec.h"
  67 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  68 #include "gromacs/mdlib/vsite.h"
  69 #include "gromacs/mdtypes/commrec.h"
  70 #include "gromacs/mdtypes/forcerec.h"
  71 #include "gromacs/mdtypes/inputrec.h"
  72 #include "gromacs/mdtypes/md_enums.h"
  73 #include "gromacs/mdtypes/mdatom.h"
  74 #include "gromacs/mdtypes/state.h"
  75 #include "gromacs/pbcutil/mshift.h"
  76 #include "gromacs/pbcutil/pbc.h"
  77 #include "gromacs/topology/mtop_util.h"
  78 #include "gromacs/topology/topsort.h"
  79 #include "gromacs/utility/cstringutil.h"
  80 #include "gromacs/utility/exceptions.h"
  81 #include "gromacs/utility/fatalerror.h"
  82 #include "gromacs/utility/gmxassert.h"
  83 #include "gromacs/utility/logger.h"
  84 #include "gromacs/utility/smalloc.h"
  85 #include "gromacs/utility/strconvert.h"
  86 #include "gromacs/utility/stringstream.h"
  87 #include "gromacs/utility/stringutil.h"
  88 #include "gromacs/utility/textwriter.h"
  89
  90 #include "domdec_constraints.h"
  91 #include "domdec_internal.h"
  92 #include "domdec_vsite.h"
  93 #include "dump.h"
  94
  95 using gmx::ArrayRef;
  96 using gmx::DDBondedChecking;
  97 using gmx::RVec;
  98
  99 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
 100 #define NITEM_DD_INIT_LOCAL_STATE 5
 101
 102 void dd_init_local_state(const gmx_domdec_t& dd, const t_state* state_global, t_state* state_local)
 103 {
 104     int buf[NITEM_DD_INIT_LOCAL_STATE];
 105
 106     if (DDMASTER(dd))
 107     {
 108         buf[0] = state_global->flags;
 109         buf[1] = state_global->ngtc;
 110         buf[2] = state_global->nnhpres;
 111         buf[3] = state_global->nhchainlength;
 112         buf[4] = state_global->dfhist ? state_global->dfhist->nlambda : 0;
 113     }
 114     dd_bcast(&dd, NITEM_DD_INIT_LOCAL_STATE * sizeof(int), buf);
 115
 116     init_gtc_state(state_local, buf[1], buf[2], buf[3]);
 117     init_dfhist_state(state_local, buf[4]);
 118     state_local->flags = buf[0];
 119 }
 120
 121 /*! \brief Check if a link is stored in \p link between charge groups \p cg_gl and \p cg_gl_j and if not so, store a link */
 122 static void check_link(t_blocka* link, int cg_gl, int cg_gl_j)
 123 {
 124     bool bFound = false;
 125     for (int k = link->index[cg_gl]; k < link->index[cg_gl + 1]; k++)
 126     {
 127         GMX_RELEASE_ASSERT(link->a, "Inconsistent NULL pointer while making charge-group links");
 128         if (link->a[k] == cg_gl_j)
 129         {
 130             bFound = TRUE;
 131         }
 132     }
 133     if (!bFound)
 134     {
 135         GMX_RELEASE_ASSERT(link->a || link->index[cg_gl + 1] + 1 > link->nalloc_a,
 136                            "Inconsistent allocation of link");
 137         /* Add this charge group link */
 138         if (link->index[cg_gl + 1] + 1 > link->nalloc_a)
 139         {
 140             link->nalloc_a = over_alloc_large(link->index[cg_gl + 1] + 1);
 141             srenew(link->a, link->nalloc_a);
 142         }
 143         link->a[link->index[cg_gl + 1]] = cg_gl_j;
 144         link->index[cg_gl + 1]++;
 145     }
 146 }
 147
 148 void makeBondedLinks(gmx_domdec_t*                              dd,
 149                      const gmx_mtop_t&                          mtop,
 150                      gmx::ArrayRef<AtomInfoWithinMoleculeBlock> atomInfoForEachMoleculeBlock)
 151 {
 152
 153     if (!dd->comm->systemInfo.filterBondedCommunication)
 154     {
 155         /* Only communicate atoms based on cut-off */
 156         dd->comm->bondedLinks = nullptr;
 157         return;
 158     }
 159
 160     t_blocka* link = nullptr;
 161
 162     /* For each atom make a list of other atoms in the system
 163      * that a linked to it via bonded interactions
 164      * which are also stored in reverse_top.
 165      */
 166
 167     reverse_ilist_t ril_intermol;
 168     if (mtop.bIntermolecularInteractions)
 169     {
 170         t_atoms atoms;
 171
 172         atoms.nr   = mtop.natoms;
 173         atoms.atom = nullptr;
 174
 175         GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
 176                            "We should have an ilist when intermolecular interactions are on");
 177
 178         ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
 179         make_reverse_ilist(
 180                 *mtop.intermolecular_ilist, &atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril_intermol);
 181     }
 182
 183     snew(link, 1);
 184     snew(link->index, mtop.natoms + 1);
 185     link->nalloc_a = 0;
 186     link->a        = nullptr;
 187
 188     link->index[0]                 = 0;
 189     int indexOfFirstAtomInMolecule = 0;
 190     int numLinkedAtoms             = 0;
 191     for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
 192     {
 193         const gmx_molblock_t& molb = mtop.molblock[mb];
 194         if (molb.nmol == 0)
 195         {
 196             continue;
 197         }
 198         const gmx_moltype_t& molt = mtop.moltype[molb.type];
 199         /* Make a reverse ilist in which the interactions are linked
 200          * to all atoms, not only the first atom as in gmx_reverse_top.
 201          * The constraints are discarded here.
 202          */
 203         ReverseTopOptions rtOptions(DDBondedChecking::ExcludeZeroLimit);
 204         reverse_ilist_t   ril;
 205         make_reverse_ilist(molt.ilist, &molt.atoms, rtOptions, AtomLinkRule::AllAtomsInBondeds, &ril);
 206
 207         AtomInfoWithinMoleculeBlock* atomInfoOfMoleculeBlock = &atomInfoForEachMoleculeBlock[mb];
 208
 209         int mol = 0;
 210         for (mol = 0; mol < (mtop.bIntermolecularInteractions ? molb.nmol : 1); mol++)
 211         {
 212             for (int a = 0; a < molt.atoms.nr; a++)
 213             {
 214                 int atomIndex              = indexOfFirstAtomInMolecule + a;
 215                 link->index[atomIndex + 1] = link->index[atomIndex];
 216                 int i                      = ril.index[a];
 217                 while (i < ril.index[a + 1])
 218                 {
 219                     int ftype = ril.il[i++];
 220                     int nral  = NRAL(ftype);
 221                     /* Skip the ifunc index */
 222                     i++;
 223                     for (int j = 0; j < nral; j++)
 224                     {
 225                         int aj = ril.il[i + j];
 226                         if (aj != a)
 227                         {
 228                             check_link(link, atomIndex, indexOfFirstAtomInMolecule + aj);
 229                         }
 230                     }
 231                     i += nral_rt(ftype);
 232                 }
 233
 234                 if (mtop.bIntermolecularInteractions)
 235                 {
 236                     int i = ril_intermol.index[atomIndex];
 237                     while (i < ril_intermol.index[atomIndex + 1])
 238                     {
 239                         int ftype = ril_intermol.il[i++];
 240                         int nral  = NRAL(ftype);
 241                         /* Skip the ifunc index */
 242                         i++;
 243                         for (int j = 0; j < nral; j++)
 244                         {
 245                             /* Here we assume we have no charge groups;
 246                              * this has been checked above.
 247                              */
 248                             int aj = ril_intermol.il[i + j];
 249                             check_link(link, atomIndex, aj);
 250                         }
 251                         i += nral_rt(ftype);
 252                     }
 253                 }
 254                 if (link->index[atomIndex + 1] - link->index[atomIndex] > 0)
 255                 {
 256                     SET_CGINFO_BOND_INTER(atomInfoOfMoleculeBlock->atomInfo[a]);
 257                     numLinkedAtoms++;
 258                 }
 259             }
 260
 261             indexOfFirstAtomInMolecule += molt.atoms.nr;
 262         }
 263         int nlink_mol = link->index[indexOfFirstAtomInMolecule]
 264                         - link->index[indexOfFirstAtomInMolecule - molt.atoms.nr];
 265
 266         if (debug)
 267         {
 268             fprintf(debug,
 269                     "molecule type '%s' %d atoms has %d atom links through bonded interac.\n",
 270                     *molt.name,
 271                     molt.atoms.nr,
 272                     nlink_mol);
 273         }
 274
 275         if (molb.nmol > mol)
 276         {
 277             /* Copy the data for the rest of the molecules in this block */
 278             link->nalloc_a += (molb.nmol - mol) * nlink_mol;
 279             srenew(link->a, link->nalloc_a);
 280             for (; mol < molb.nmol; mol++)
 281             {
 282                 for (int a = 0; a < molt.atoms.nr; a++)
 283                 {
 284                     int atomIndex              = indexOfFirstAtomInMolecule + a;
 285                     link->index[atomIndex + 1] = link->index[atomIndex + 1 - molt.atoms.nr] + nlink_mol;
 286                     for (int j = link->index[atomIndex]; j < link->index[atomIndex + 1]; j++)
 287                     {
 288                         link->a[j] = link->a[j - nlink_mol] + molt.atoms.nr;
 289                     }
 290                     if (link->index[atomIndex + 1] - link->index[atomIndex] > 0
 291                         && atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock
 292                                    < gmx::ssize(atomInfoOfMoleculeBlock->atomInfo))
 293                     {
 294                         SET_CGINFO_BOND_INTER(
 295                                 atomInfoOfMoleculeBlock
 296                                         ->atomInfo[atomIndex - atomInfoOfMoleculeBlock->indexOfFirstAtomInMoleculeBlock]);
 297                         numLinkedAtoms++;
 298                     }
 299                 }
 300                 indexOfFirstAtomInMolecule += molt.atoms.nr;
 301             }
 302         }
 303     }
 304
 305     if (debug)
 306     {
 307         fprintf(debug, "Of the %d atoms %d are linked via bonded interactions\n", mtop.natoms, numLinkedAtoms);
 308     }
 309
 310     dd->comm->bondedLinks = link;
 311 }
 312
 313 typedef struct
 314 {
 315     real r2;
 316     int  ftype;
 317     int  a1;
 318     int  a2;
 319 } bonded_distance_t;
 320
 321 /*! \brief Compare distance^2 \p r2 against the distance in \p bd and if larger store it along with \p ftype and atom indices \p a1 and \p a2 */
 322 static void update_max_bonded_distance(real r2, int ftype, int a1, int a2, bonded_distance_t* bd)
 323 {
 324     if (r2 > bd->r2)
 325     {
 326         bd->r2    = r2;
 327         bd->ftype = ftype;
 328         bd->a1    = a1;
 329         bd->a2    = a2;
 330     }
 331 }
 332
 333 /*! \brief Set the distance, function type and atom indices for the longest distance between atoms of molecule type \p molt for two-body and multi-body bonded interactions */
 334 static void bonded_cg_distance_mol(const gmx_moltype_t*   molt,
 335                                    const DDBondedChecking ddBondedChecking,
 336                                    gmx_bool               bExcl,
 337                                    ArrayRef<const RVec>   x,
 338                                    bonded_distance_t*     bd_2b,
 339                                    bonded_distance_t*     bd_mb)
 340 {
 341     const ReverseTopOptions rtOptions(ddBondedChecking);
 342
 343     for (int ftype = 0; ftype < F_NRE; ftype++)
 344     {
 345         if (dd_check_ftype(ftype, rtOptions))
 346         {
 347             const auto& il   = molt->ilist[ftype];
 348             int         nral = NRAL(ftype);
 349             if (nral > 1)
 350             {
 351                 for (int i = 0; i < il.size(); i += 1 + nral)
 352                 {
 353                     for (int ai = 0; ai < nral; ai++)
 354                     {
 355                         int atomI = il.iatoms[i + 1 + ai];
 356                         for (int aj = ai + 1; aj < nral; aj++)
 357                         {
 358                             int atomJ = il.iatoms[i + 1 + aj];
 359                             if (atomI != atomJ)
 360                             {
 361                                 real rij2 = distance2(x[atomI], x[atomJ]);
 362
 363                                 update_max_bonded_distance(
 364                                         rij2, ftype, atomI, atomJ, (nral == 2) ? bd_2b : bd_mb);
 365                             }
 366                         }
 367                     }
 368                 }
 369             }
 370         }
 371     }
 372     if (bExcl)
 373     {
 374         const auto& excls = molt->excls;
 375         for (gmx::index ai = 0; ai < excls.ssize(); ai++)
 376         {
 377             for (const int aj : excls[ai])
 378             {
 379                 if (ai != aj)
 380                 {
 381                     real rij2 = distance2(x[ai], x[aj]);
 382
 383                     /* There is no function type for exclusions, use -1 */
 384                     update_max_bonded_distance(rij2, -1, ai, aj, bd_2b);
 385                 }
 386             }
 387         }
 388     }
 389 }
 390
 391 /*! \brief Set the distance, function type and atom indices for the longest atom distance involved in intermolecular interactions for two-body and multi-body bonded interactions */
 392 static void bonded_distance_intermol(const InteractionLists& ilists_intermol,
 393                                      const DDBondedChecking  ddBondedChecking,
 394                                      ArrayRef<const RVec>    x,
 395                                      PbcType                 pbcType,
 396                                      const matrix            box,
 397                                      bonded_distance_t*      bd_2b,
 398                                      bonded_distance_t*      bd_mb)
 399 {
 400     t_pbc pbc;
 401
 402     set_pbc(&pbc, pbcType, box);
 403
 404     const ReverseTopOptions rtOptions(ddBondedChecking);
 405
 406     for (int ftype = 0; ftype < F_NRE; ftype++)
 407     {
 408         if (dd_check_ftype(ftype, rtOptions))
 409         {
 410             const auto& il   = ilists_intermol[ftype];
 411             int         nral = NRAL(ftype);
 412
 413             /* No nral>1 check here, since intermol interactions always
 414              * have nral>=2 (and the code is also correct for nral=1).
 415              */
 416             for (int i = 0; i < il.size(); i += 1 + nral)
 417             {
 418                 for (int ai = 0; ai < nral; ai++)
 419                 {
 420                     int atom_i = il.iatoms[i + 1 + ai];
 421
 422                     for (int aj = ai + 1; aj < nral; aj++)
 423                     {
 424                         rvec dx;
 425
 426                         int atom_j = il.iatoms[i + 1 + aj];
 427
 428                         pbc_dx(&pbc, x[atom_i], x[atom_j], dx);
 429
 430                         const real rij2 = norm2(dx);
 431
 432                         update_max_bonded_distance(rij2, ftype, atom_i, atom_j, (nral == 2) ? bd_2b : bd_mb);
 433                     }
 434                 }
 435             }
 436         }
 437     }
 438 }
 439
 440 //! Returns whether \p molt has at least one virtual site
 441 static bool moltypeHasVsite(const gmx_moltype_t& molt)
 442 {
 443     bool hasVsite = false;
 444     for (int i = 0; i < F_NRE; i++)
 445     {
 446         if ((interaction_function[i].flags & IF_VSITE) && !molt.ilist[i].empty())
 447         {
 448             hasVsite = true;
 449         }
 450     }
 451
 452     return hasVsite;
 453 }
 454
 455 //! Returns coordinates not broken over PBC for a molecule
 456 static void getWholeMoleculeCoordinates(const gmx_moltype_t*  molt,
 457                                         const gmx_ffparams_t* ffparams,
 458                                         PbcType               pbcType,
 459                                         t_graph*              graph,
 460                                         const matrix          box,
 461                                         ArrayRef<const RVec>  x,
 462                                         ArrayRef<RVec>        xs)
 463 {
 464     if (pbcType != PbcType::No)
 465     {
 466         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(x.data()));
 467
 468         shift_x(graph, box, as_rvec_array(x.data()), as_rvec_array(xs.data()));
 469         /* By doing an extra mk_mshift the molecules that are broken
 470          * because they were e.g. imported from another software
 471          * will be made whole again. Such are the healing powers
 472          * of GROMACS.
 473          */
 474         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(xs.data()));
 475     }
 476     else
 477     {
 478         /* We copy the coordinates so the original coordinates remain
 479          * unchanged, just to be 100% sure that we do not affect
 480          * binary reproducibility of simulations.
 481          */
 482         for (int i = 0; i < molt->atoms.nr; i++)
 483         {
 484             copy_rvec(x[i], xs[i]);
 485         }
 486     }
 487
 488     if (moltypeHasVsite(*molt))
 489     {
 490         gmx::constructVirtualSites(xs, ffparams->iparams, molt->ilist);
 491     }
 492 }
 493
 494 void dd_bonded_cg_distance(const gmx::MDLogger&   mdlog,
 495                            const gmx_mtop_t&      mtop,
 496                            const t_inputrec&      inputrec,
 497                            ArrayRef<const RVec>   x,
 498                            const matrix           box,
 499                            const DDBondedChecking ddBondedChecking,
 500                            real*                  r_2b,
 501                            real*                  r_mb)
 502 {
 503     bonded_distance_t bd_2b = { 0, -1, -1, -1 };
 504     bonded_distance_t bd_mb = { 0, -1, -1, -1 };
 505
 506     bool bExclRequired = inputrecExclForces(&inputrec);
 507
 508     *r_2b         = 0;
 509     *r_mb         = 0;
 510     int at_offset = 0;
 511     for (const gmx_molblock_t& molb : mtop.molblock)
 512     {
 513         const gmx_moltype_t& molt = mtop.moltype[molb.type];
 514         if (molt.atoms.nr == 1 || molb.nmol == 0)
 515         {
 516             at_offset += molb.nmol * molt.atoms.nr;
 517         }
 518         else
 519         {
 520             t_graph graph;
 521             if (inputrec.pbcType != PbcType::No)
 522             {
 523                 graph = mk_graph_moltype(molt);
 524             }
 525
 526             std::vector<RVec> xs(molt.atoms.nr);
 527             for (int mol = 0; mol < molb.nmol; mol++)
 528             {
 529                 getWholeMoleculeCoordinates(&molt,
 530                                             &mtop.ffparams,
 531                                             inputrec.pbcType,
 532                                             &graph,
 533                                             box,
 534                                             x.subArray(at_offset, molt.atoms.nr),
 535                                             xs);
 536
 537                 bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 };
 538                 bonded_distance_t bd_mol_mb = { 0, -1, -1, -1 };
 539
 540                 bonded_cg_distance_mol(&molt, ddBondedChecking, bExclRequired, xs, &bd_mol_2b, &bd_mol_mb);
 541
 542                 /* Process the mol data adding the atom index offset */
 543                 update_max_bonded_distance(bd_mol_2b.r2,
 544                                            bd_mol_2b.ftype,
 545                                            at_offset + bd_mol_2b.a1,
 546                                            at_offset + bd_mol_2b.a2,
 547                                            &bd_2b);
 548                 update_max_bonded_distance(bd_mol_mb.r2,
 549                                            bd_mol_mb.ftype,
 550                                            at_offset + bd_mol_mb.a1,
 551                                            at_offset + bd_mol_mb.a2,
 552                                            &bd_mb);
 553
 554                 at_offset += molt.atoms.nr;
 555             }
 556         }
 557     }
 558
 559     if (mtop.bIntermolecularInteractions)
 560     {
 561         GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
 562                            "We should have an ilist when intermolecular interactions are on");
 563
 564         bonded_distance_intermol(
 565                 *mtop.intermolecular_ilist, ddBondedChecking, x, inputrec.pbcType, box, &bd_2b, &bd_mb);
 566     }
 567
 568     *r_2b = sqrt(bd_2b.r2);
 569     *r_mb = sqrt(bd_mb.r2);
 570
 571     if (*r_2b > 0 || *r_mb > 0)
 572     {
 573         GMX_LOG(mdlog.info).appendText("Initial maximum distances in bonded interactions:");
 574         if (*r_2b > 0)
 575         {
 576             GMX_LOG(mdlog.info)
 577                     .appendTextFormatted(
 578                             "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d",
 579                             *r_2b,
 580                             (bd_2b.ftype >= 0) ? interaction_function[bd_2b.ftype].longname : "Exclusion",
 581                             bd_2b.a1 + 1,
 582                             bd_2b.a2 + 1);
 583         }
 584         if (*r_mb > 0)
 585         {
 586             GMX_LOG(mdlog.info)
 587                     .appendTextFormatted(
 588                             "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d",
 589                             *r_mb,
 590                             interaction_function[bd_mb.ftype].longname,
 591                             bd_mb.a1 + 1,
 592                             bd_mb.a2 + 1);
 593         }
 594     }
 595 }