src/gromacs/domdec/localtopologychecker.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  *
  38  * \brief This file defines functions used by the domdec module
  39  * while managing the construction, use and error checking for
  40  * topologies local to a DD rank.
  41  *
  42  * \ingroup module_domdec
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include <string>
  48 #include <vector>
  49
  50 #include "gromacs/domdec/localtopologychecker.h"
  51
  52 #include "gromacs/domdec/domdec_internal.h"
  53 #include "gromacs/domdec/reversetopology.h"
  54 #include "gromacs/gmxlib/network.h"
  55 #include "gromacs/mdtypes/commrec.h"
  56 #include "gromacs/topology/idef.h"
  57 #include "gromacs/topology/ifunc.h"
  58 #include "gromacs/topology/mtop_util.h"
  59 #include "gromacs/utility/arrayref.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/logger.h"
  62 #include "gromacs/utility/stringstream.h"
  63 #include "gromacs/utility/textwriter.h"
  64
  65 #include "dump.h"
  66
  67 using gmx::ArrayRef;
  68 using gmx::RVec;
  69
  70 /*! \brief Checks whether interactions have been assigned for one function type
  71  *
  72  * Loops over a list of interactions in the local topology of one function type
  73  * and flags each of the interactions as assigned in the global \p isAssigned list.
  74  * Exits with an inconsistency error when an interaction is assigned more than once.
  75  */
  76 static void flagInteractionsForType(const int              ftype,
  77                                     const InteractionList& il,
  78                                     const reverse_ilist_t& ril,
  79                                     const gmx::Range<int>& atomRange,
  80                                     const int              numAtomsPerMolecule,
  81                                     ArrayRef<const int>    globalAtomIndices,
  82                                     ArrayRef<int>          isAssigned)
  83 {
  84     const int nril_mol = ril.index[numAtomsPerMolecule];
  85     const int nral     = NRAL(ftype);
  86
  87     for (int i = 0; i < il.size(); i += 1 + nral)
  88     {
  89         // ia[0] is the interaction type, ia[1, ...] the atom indices
  90         const int* ia = il.iatoms.data() + i;
  91         // Extract the global atom index of the first atom in this interaction
  92         const int a0 = globalAtomIndices[ia[1]];
  93         /* Check if this interaction is in
  94          * the currently checked molblock.
  95          */
  96         if (atomRange.isInRange(a0))
  97         {
  98             // The molecule index in the list of this molecule type
  99             const int moleculeIndex = (a0 - atomRange.begin()) / numAtomsPerMolecule;
 100             const int atomOffset = (a0 - atomRange.begin()) - moleculeIndex * numAtomsPerMolecule;
 101             const int globalAtomStartInMolecule = atomRange.begin() + moleculeIndex * numAtomsPerMolecule;
 102             int       j_mol                     = ril.index[atomOffset];
 103             bool found                          = false;
 104             while (j_mol < ril.index[atomOffset + 1] && !found)
 105             {
 106                 const int j       = moleculeIndex * nril_mol + j_mol;
 107                 const int ftype_j = ril.il[j_mol];
 108                 /* Here we need to check if this interaction has
 109                  * not already been assigned, since we could have
 110                  * multiply defined interactions.
 111                  */
 112                 if (ftype == ftype_j && ia[0] == ril.il[j_mol + 1] && isAssigned[j] == 0)
 113                 {
 114                     /* Check the atoms */
 115                     found = true;
 116                     for (int a = 0; a < nral; a++)
 117                     {
 118                         if (globalAtomIndices[ia[1 + a]]
 119                             != globalAtomStartInMolecule + ril.il[j_mol + 2 + a])
 120                         {
 121                             found = false;
 122                         }
 123                     }
 124                     if (found)
 125                     {
 126                         isAssigned[j] = 1;
 127                     }
 128                 }
 129                 j_mol += 2 + nral_rt(ftype_j);
 130             }
 131             if (!found)
 132             {
 133                 gmx_incons("Some interactions seem to be assigned multiple times");
 134             }
 135         }
 136     }
 137 }
 138
 139 /*! \brief Help print error output when interactions are missing in a molblock
 140  *
 141  * \note This function needs to be called on all ranks (contains a global summation)
 142  */
 143 static std::string printMissingInteractionsMolblock(t_commrec*               cr,
 144                                                     const gmx_reverse_top_t& rt,
 145                                                     const char*              moltypename,
 146                                                     const reverse_ilist_t&   ril,
 147                                                     const gmx::Range<int>&   atomRange,
 148                                                     const int                numAtomsPerMolecule,
 149                                                     const int                numMolecules,
 150                                                     const InteractionDefinitions& idef)
 151 {
 152     const int               nril_mol = ril.index[numAtomsPerMolecule];
 153     std::vector<int>        isAssigned(numMolecules * nril_mol, 0);
 154     gmx::StringOutputStream stream;
 155     gmx::TextWriter         log(&stream);
 156
 157     for (int ftype = 0; ftype < F_NRE; ftype++)
 158     {
 159         if (dd_check_ftype(ftype, rt.options()))
 160         {
 161             flagInteractionsForType(
 162                     ftype, idef.il[ftype], ril, atomRange, numAtomsPerMolecule, cr->dd->globalAtomIndices, isAssigned);
 163         }
 164     }
 165
 166     gmx_sumi(isAssigned.size(), isAssigned.data(), cr);
 167
 168     const int numMissingToPrint = 10;
 169     int       i                 = 0;
 170     for (int mol = 0; mol < numMolecules; mol++)
 171     {
 172         int j_mol = 0;
 173         while (j_mol < nril_mol)
 174         {
 175             int ftype = ril.il[j_mol];
 176             int nral  = NRAL(ftype);
 177             int j     = mol * nril_mol + j_mol;
 178             if (isAssigned[j] == 0 && !(interaction_function[ftype].flags & IF_VSITE))
 179             {
 180                 if (DDMASTER(cr->dd))
 181                 {
 182                     if (i == 0)
 183                     {
 184                         log.writeLineFormatted("Molecule type '%s'", moltypename);
 185                         log.writeLineFormatted(
 186                                 "the first %d missing interactions, except for exclusions:",
 187                                 numMissingToPrint);
 188                     }
 189                     log.writeStringFormatted("%20s atoms", interaction_function[ftype].longname);
 190                     int a = 0;
 191                     for (; a < nral; a++)
 192                     {
 193                         log.writeStringFormatted("%5d", ril.il[j_mol + 2 + a] + 1);
 194                     }
 195                     while (a < 4)
 196                     {
 197                         log.writeString("     ");
 198                         a++;
 199                     }
 200                     log.writeString(" global");
 201                     for (int a = 0; a < nral; a++)
 202                     {
 203                         log.writeStringFormatted("%6d",
 204                                                  atomRange.begin() + mol * numAtomsPerMolecule
 205                                                          + ril.il[j_mol + 2 + a] + 1);
 206                     }
 207                     log.ensureLineBreak();
 208                 }
 209                 i++;
 210                 if (i >= numMissingToPrint)
 211                 {
 212                     break;
 213                 }
 214             }
 215             j_mol += 2 + nral_rt(ftype);
 216         }
 217     }
 218
 219     return stream.toString();
 220 }
 221
 222 /*! \brief Help print error output when interactions are missing */
 223 static void printMissingInteractionsAtoms(const gmx::MDLogger&          mdlog,
 224                                           t_commrec*                    cr,
 225                                           const gmx_mtop_t&             mtop,
 226                                           const InteractionDefinitions& idef)
 227 {
 228     const gmx_reverse_top_t& rt = *cr->dd->reverse_top;
 229
 230     /* Print the atoms in the missing interactions per molblock */
 231     int a_end = 0;
 232     for (const gmx_molblock_t& molb : mtop.molblock)
 233     {
 234         const gmx_moltype_t& moltype = mtop.moltype[molb.type];
 235         const int            a_start = a_end;
 236         a_end                        = a_start + molb.nmol * moltype.atoms.nr;
 237         const gmx::Range<int> atomRange(a_start, a_end);
 238
 239         auto warning = printMissingInteractionsMolblock(
 240                 cr,
 241                 rt,
 242                 *(moltype.name),
 243                 cr->dd->reverse_top->interactionListForMoleculeType(molb.type),
 244                 atomRange,
 245                 moltype.atoms.nr,
 246                 molb.nmol,
 247                 idef);
 248
 249         GMX_LOG(mdlog.warning).appendText(warning);
 250     }
 251 }
 252
 253 /*! \brief Print error output when interactions are missing */
 254 [[noreturn]] static void dd_print_missing_interactions(const gmx::MDLogger& mdlog,
 255                                                        t_commrec*           cr,
 256                                                        const int numBondedInteractionsOverAllDomains,
 257                                                        const int expectedNumGlobalBondedInteractions,
 258                                                        const gmx_mtop_t&     top_global,
 259                                                        const gmx_localtop_t* top_local,
 260                                                        ArrayRef<const RVec>  x,
 261                                                        const matrix          box)
 262 {
 263     int           cl[F_NRE];
 264     gmx_domdec_t* dd = cr->dd;
 265
 266     GMX_LOG(mdlog.warning)
 267             .appendText(
 268                     "Not all bonded interactions have been properly assigned to the domain "
 269                     "decomposition cells");
 270
 271     const int ndiff_tot = numBondedInteractionsOverAllDomains - expectedNumGlobalBondedInteractions;
 272
 273     for (int ftype = 0; ftype < F_NRE; ftype++)
 274     {
 275         const int nral = NRAL(ftype);
 276         cl[ftype]      = top_local->idef.il[ftype].size() / (1 + nral);
 277     }
 278
 279     gmx_sumi(F_NRE, cl, cr);
 280
 281     if (DDMASTER(dd))
 282     {
 283         GMX_LOG(mdlog.warning).appendText("A list of missing interactions:");
 284         int rest_global = expectedNumGlobalBondedInteractions;
 285         int rest        = numBondedInteractionsOverAllDomains;
 286         for (int ftype = 0; ftype < F_NRE; ftype++)
 287         {
 288             /* In the reverse and local top all constraints are merged
 289              * into F_CONSTR. So in the if statement we skip F_CONSTRNC
 290              * and add these constraints when doing F_CONSTR.
 291              */
 292             if (dd_check_ftype(ftype, dd->reverse_top->options()) && ftype != F_CONSTRNC)
 293             {
 294                 int n = gmx_mtop_ftype_count(top_global, ftype);
 295                 if (ftype == F_CONSTR)
 296                 {
 297                     n += gmx_mtop_ftype_count(top_global, F_CONSTRNC);
 298                 }
 299                 int ndiff = cl[ftype] - n;
 300                 if (ndiff != 0)
 301                 {
 302                     GMX_LOG(mdlog.warning)
 303                             .appendTextFormatted("%20s of %6d missing %6d",
 304                                                  interaction_function[ftype].longname,
 305                                                  n,
 306                                                  -ndiff);
 307                 }
 308                 rest_global -= n;
 309                 rest -= cl[ftype];
 310             }
 311         }
 312
 313         int ndiff = rest - rest_global;
 314         if (ndiff != 0)
 315         {
 316             GMX_LOG(mdlog.warning).appendTextFormatted("%20s of %6d missing %6d", "exclusions", rest_global, -ndiff);
 317         }
 318     }
 319
 320     printMissingInteractionsAtoms(mdlog, cr, top_global, top_local->idef);
 321     write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr, -1, as_rvec_array(x.data()), box);
 322
 323     std::string errorMessage;
 324
 325     if (ndiff_tot > 0)
 326     {
 327         errorMessage =
 328                 "One or more interactions were assigned to multiple domains of the domain "
 329                 "decompostion. Please report this bug.";
 330     }
 331     else
 332     {
 333         errorMessage = gmx::formatString(
 334                 "%d of the %d bonded interactions could not be calculated because some atoms "
 335                 "involved moved further apart than the multi-body cut-off distance (%g nm) or the "
 336                 "two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds "
 337                 "also see option -ddcheck",
 338                 -ndiff_tot,
 339                 expectedNumGlobalBondedInteractions,
 340                 dd_cutoff_multibody(dd),
 341                 dd_cutoff_twobody(dd));
 342     }
 343     gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), "%s", errorMessage.c_str());
 344 }
 345
 346 namespace gmx
 347 {
 348
 349 /*! \brief Compute the total bonded interaction count
 350  *
 351  * \param[in] mtop              The global system topology
 352  * \param[in] useUpdateGroups   Whether update groups are in use
 353  *
 354  * When using domain decomposition without update groups,
 355  * constraint-type interactions can be split across domains, and so we
 356  * do not consider them in this correctness check. Otherwise, we
 357  * include them.
 358  */
 359 static int computeExpectedNumGlobalBondedInteractions(const gmx_mtop_t& mtop, const bool useUpdateGroups)
 360 {
 361     int expectedNumGlobalBondedInteractions = gmx_mtop_interaction_count(mtop, IF_BOND);
 362     if (useUpdateGroups)
 363     {
 364         expectedNumGlobalBondedInteractions += gmx_mtop_interaction_count(mtop, IF_CONSTRAINT);
 365     }
 366     return expectedNumGlobalBondedInteractions;
 367 }
 368
 369 LocalTopologyChecker::LocalTopologyChecker(const gmx_mtop_t& mtop, const bool useUpdateGroups) :
 370     expectedNumGlobalBondedInteractions(computeExpectedNumGlobalBondedInteractions(mtop, useUpdateGroups))
 371 {
 372 }
 373
 374 } // namespace gmx
 375
 376 void scheduleCheckOfLocalTopology(gmx_domdec_t* dd, const int numBondedInteractionsToReduce)
 377 {
 378     dd->localTopologyChecker->numBondedInteractionsToReduce = numBondedInteractionsToReduce;
 379     // Note that it's possible for this to still be true from the last
 380     // time it was set, e.g. if repartitioning was triggered before
 381     // global communication that would have acted on the true
 382     // value. This could happen for example when replica exchange took
 383     // place soon after a partition.
 384     dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = true;
 385     // Clear the old global value, which is now invalid
 386     dd->localTopologyChecker->numBondedInteractionsOverAllDomains.reset();
 387 }
 388
 389 bool shouldCheckNumberOfBondedInteractions(const gmx_domdec_t& dd)
 390 {
 391     return dd.localTopologyChecker->shouldCheckNumberOfBondedInteractions;
 392 }
 393
 394 int numBondedInteractions(const gmx_domdec_t& dd)
 395 {
 396     return dd.localTopologyChecker->numBondedInteractionsToReduce;
 397 }
 398
 399 void setNumberOfBondedInteractionsOverAllDomains(gmx_domdec_t* dd, int newValue)
 400 {
 401     GMX_RELEASE_ASSERT(!dd->localTopologyChecker->numBondedInteractionsOverAllDomains.has_value(),
 402                        "Cannot set number of bonded interactions because it is already set");
 403     dd->localTopologyChecker->numBondedInteractionsOverAllDomains.emplace(newValue);
 404 }
 405
 406 void checkNumberOfBondedInteractions(const gmx::MDLogger&  mdlog,
 407                                      t_commrec*            cr,
 408                                      const gmx_mtop_t&     top_global,
 409                                      const gmx_localtop_t* top_local,
 410                                      ArrayRef<const RVec>  x,
 411                                      const matrix          box)
 412 {
 413     GMX_RELEASE_ASSERT(
 414             DOMAINDECOMP(cr),
 415             "No need to check number of bonded interactions when not using domain decomposition");
 416     if (cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions)
 417     {
 418         GMX_RELEASE_ASSERT(cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.has_value(),
 419                            "The check for the total number of bonded interactions requires the "
 420                            "value to have been reduced across all domains");
 421         if (cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value()
 422             != cr->dd->localTopologyChecker->expectedNumGlobalBondedInteractions)
 423         {
 424             dd_print_missing_interactions(
 425                     mdlog,
 426                     cr,
 427                     cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value(),
 428                     cr->dd->localTopologyChecker->expectedNumGlobalBondedInteractions,
 429                     top_global,
 430                     top_local,
 431                     x,
 432                     box); // Does not return
 433         }
 434         // Now that the value is set and the check complete, future
 435         // global communication should not compute the value until
 436         // after the next partitioning.
 437         cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = false;
 438     }
 439 }