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 "gromacs/domdec/localtopologychecker.h"
  48
  49 #include <optional>
  50 #include <string>
  51 #include <vector>
  52
  53 #include "gromacs/domdec/domdec_internal.h"
  54 #include "gromacs/domdec/reversetopology.h"
  55 #include "gromacs/gmxlib/network.h"
  56 #include "gromacs/mdtypes/commrec.h"
  57 #include "gromacs/topology/idef.h"
  58 #include "gromacs/topology/ifunc.h"
  59 #include "gromacs/topology/mtop_util.h"
  60 #include "gromacs/utility/arrayref.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/logger.h"
  63 #include "gromacs/utility/stringstream.h"
  64 #include "gromacs/utility/textwriter.h"
  65
  66 #include "dump.h"
  67
  68 namespace gmx
  69 {
  70
  71 /*! \brief Checks whether interactions have been assigned for one function type
  72  *
  73  * Loops over a list of interactions in the local topology of one function type
  74  * and flags each of the interactions as assigned in the global \p isAssigned list.
  75  * Exits with an inconsistency error when an interaction is assigned more than once.
  76  */
  77 static void flagInteractionsForType(const int              ftype,
  78                                     const InteractionList& il,
  79                                     const reverse_ilist_t& ril,
  80                                     const Range<int>&      atomRange,
  81                                     const int              numAtomsPerMolecule,
  82                                     ArrayRef<const int>    globalAtomIndices,
  83                                     ArrayRef<int>          isAssigned)
  84 {
  85     const int nril_mol = ril.index[numAtomsPerMolecule];
  86     const int nral     = NRAL(ftype);
  87
  88     for (int i = 0; i < il.size(); i += 1 + nral)
  89     {
  90         // ia[0] is the interaction type, ia[1, ...] the atom indices
  91         const int* ia = il.iatoms.data() + i;
  92         // Extract the global atom index of the first atom in this interaction
  93         const int a0 = globalAtomIndices[ia[1]];
  94         /* Check if this interaction is in
  95          * the currently checked molblock.
  96          */
  97         if (atomRange.isInRange(a0))
  98         {
  99             // The molecule index in the list of this molecule type
 100             const int moleculeIndex = (a0 - atomRange.begin()) / numAtomsPerMolecule;
 101             const int atomOffset = (a0 - atomRange.begin()) - moleculeIndex * numAtomsPerMolecule;
 102             const int globalAtomStartInMolecule = atomRange.begin() + moleculeIndex * numAtomsPerMolecule;
 103             int       j_mol                     = ril.index[atomOffset];
 104             bool found                          = false;
 105             while (j_mol < ril.index[atomOffset + 1] && !found)
 106             {
 107                 const int j       = moleculeIndex * nril_mol + j_mol;
 108                 const int ftype_j = ril.il[j_mol];
 109                 /* Here we need to check if this interaction has
 110                  * not already been assigned, since we could have
 111                  * multiply defined interactions.
 112                  */
 113                 if (ftype == ftype_j && ia[0] == ril.il[j_mol + 1] && isAssigned[j] == 0)
 114                 {
 115                     /* Check the atoms */
 116                     found = true;
 117                     for (int a = 0; a < nral; a++)
 118                     {
 119                         if (globalAtomIndices[ia[1 + a]]
 120                             != globalAtomStartInMolecule + ril.il[j_mol + 2 + a])
 121                         {
 122                             found = false;
 123                         }
 124                     }
 125                     if (found)
 126                     {
 127                         isAssigned[j] = 1;
 128                     }
 129                 }
 130                 j_mol += 2 + nral_rt(ftype_j);
 131             }
 132             if (!found)
 133             {
 134                 gmx_incons("Some interactions seem to be assigned multiple times");
 135             }
 136         }
 137     }
 138 }
 139
 140 /*! \brief Help print error output when interactions are missing in a molblock
 141  *
 142  * \note This function needs to be called on all ranks (contains a global summation)
 143  */
 144 static std::string printMissingInteractionsMolblock(const t_commrec*         cr,
 145                                                     const gmx_reverse_top_t& rt,
 146                                                     const char*              moltypename,
 147                                                     const reverse_ilist_t&   ril,
 148                                                     const Range<int>&        atomRange,
 149                                                     const int                numAtomsPerMolecule,
 150                                                     const int                numMolecules,
 151                                                     const InteractionDefinitions& idef)
 152 {
 153     const int          nril_mol = ril.index[numAtomsPerMolecule];
 154     std::vector<int>   isAssigned(numMolecules * nril_mol, 0);
 155     StringOutputStream stream;
 156     TextWriter         log(&stream);
 157
 158     for (int ftype = 0; ftype < F_NRE; ftype++)
 159     {
 160         if (dd_check_ftype(ftype, rt.options()))
 161         {
 162             flagInteractionsForType(
 163                     ftype, idef.il[ftype], ril, atomRange, numAtomsPerMolecule, cr->dd->globalAtomIndices, isAssigned);
 164         }
 165     }
 166
 167     gmx_sumi(isAssigned.size(), isAssigned.data(), cr);
 168
 169     const int numMissingToPrint = 10;
 170     int       i                 = 0;
 171     for (int mol = 0; mol < numMolecules; mol++)
 172     {
 173         int j_mol = 0;
 174         while (j_mol < nril_mol)
 175         {
 176             int ftype = ril.il[j_mol];
 177             int nral  = NRAL(ftype);
 178             int j     = mol * nril_mol + j_mol;
 179             if (isAssigned[j] == 0 && !(interaction_function[ftype].flags & IF_VSITE))
 180             {
 181                 if (DDMASTER(cr->dd))
 182                 {
 183                     if (i == 0)
 184                     {
 185                         log.writeLineFormatted("Molecule type '%s'", moltypename);
 186                         log.writeLineFormatted(
 187                                 "the first %d missing interactions, except for exclusions:",
 188                                 numMissingToPrint);
 189                     }
 190                     log.writeStringFormatted("%20s atoms", interaction_function[ftype].longname);
 191                     int a = 0;
 192                     for (; a < nral; a++)
 193                     {
 194                         log.writeStringFormatted("%5d", ril.il[j_mol + 2 + a] + 1);
 195                     }
 196                     while (a < 4)
 197                     {
 198                         log.writeString("     ");
 199                         a++;
 200                     }
 201                     log.writeString(" global");
 202                     for (int a = 0; a < nral; a++)
 203                     {
 204                         log.writeStringFormatted("%6d",
 205                                                  atomRange.begin() + mol * numAtomsPerMolecule
 206                                                          + ril.il[j_mol + 2 + a] + 1);
 207                     }
 208                     log.ensureLineBreak();
 209                 }
 210                 i++;
 211                 if (i >= numMissingToPrint)
 212                 {
 213                     break;
 214                 }
 215             }
 216             j_mol += 2 + nral_rt(ftype);
 217         }
 218     }
 219
 220     return stream.toString();
 221 }
 222
 223 /*! \brief Help print error output when interactions are missing */
 224 static void printMissingInteractionsAtoms(const MDLogger&               mdlog,
 225                                           const t_commrec*              cr,
 226                                           const gmx_mtop_t&             mtop,
 227                                           const InteractionDefinitions& idef)
 228 {
 229     const gmx_reverse_top_t& rt = *cr->dd->reverse_top;
 230
 231     /* Print the atoms in the missing interactions per molblock */
 232     int a_end = 0;
 233     for (const gmx_molblock_t& molb : mtop.molblock)
 234     {
 235         const gmx_moltype_t& moltype = mtop.moltype[molb.type];
 236         const int            a_start = a_end;
 237         a_end                        = a_start + molb.nmol * moltype.atoms.nr;
 238         const Range<int> atomRange(a_start, a_end);
 239
 240         auto warning = printMissingInteractionsMolblock(
 241                 cr,
 242                 rt,
 243                 *(moltype.name),
 244                 cr->dd->reverse_top->interactionListForMoleculeType(molb.type),
 245                 atomRange,
 246                 moltype.atoms.nr,
 247                 molb.nmol,
 248                 idef);
 249
 250         GMX_LOG(mdlog.warning).appendText(warning);
 251     }
 252 }
 253
 254 /*! \brief Print error output when interactions are missing */
 255 [[noreturn]] static void dd_print_missing_interactions(const MDLogger&  mdlog,
 256                                                        const t_commrec* cr,
 257                                                        const int numBondedInteractionsOverAllDomains,
 258                                                        const int expectedNumGlobalBondedInteractions,
 259                                                        const gmx_mtop_t&     top_global,
 260                                                        const gmx_localtop_t* top_local,
 261                                                        ArrayRef<const RVec>  x,
 262                                                        const matrix          box)
 263 {
 264     int           cl[F_NRE];
 265     gmx_domdec_t* dd = cr->dd;
 266
 267     GMX_LOG(mdlog.warning)
 268             .appendText(
 269                     "Not all bonded interactions have been properly assigned to the domain "
 270                     "decomposition cells");
 271
 272     const int ndiff_tot = numBondedInteractionsOverAllDomains - expectedNumGlobalBondedInteractions;
 273
 274     for (int ftype = 0; ftype < F_NRE; ftype++)
 275     {
 276         const int nral = NRAL(ftype);
 277         cl[ftype]      = top_local->idef.il[ftype].size() / (1 + nral);
 278     }
 279
 280     gmx_sumi(F_NRE, cl, cr);
 281
 282     if (DDMASTER(dd))
 283     {
 284         GMX_LOG(mdlog.warning).appendText("A list of missing interactions:");
 285         int rest_global = expectedNumGlobalBondedInteractions;
 286         int rest        = numBondedInteractionsOverAllDomains;
 287         for (int ftype = 0; ftype < F_NRE; ftype++)
 288         {
 289             /* In the reverse and local top all constraints are merged
 290              * into F_CONSTR. So in the if statement we skip F_CONSTRNC
 291              * and add these constraints when doing F_CONSTR.
 292              */
 293             if (dd_check_ftype(ftype, dd->reverse_top->options()) && ftype != F_CONSTRNC)
 294             {
 295                 int n = gmx_mtop_ftype_count(top_global, ftype);
 296                 if (ftype == F_CONSTR)
 297                 {
 298                     n += gmx_mtop_ftype_count(top_global, F_CONSTRNC);
 299                 }
 300                 int ndiff = cl[ftype] - n;
 301                 if (ndiff != 0)
 302                 {
 303                     GMX_LOG(mdlog.warning)
 304                             .appendTextFormatted("%20s of %6d missing %6d",
 305                                                  interaction_function[ftype].longname,
 306                                                  n,
 307                                                  -ndiff);
 308                 }
 309                 rest_global -= n;
 310                 rest -= cl[ftype];
 311             }
 312         }
 313
 314         int ndiff = rest - rest_global;
 315         if (ndiff != 0)
 316         {
 317             GMX_LOG(mdlog.warning).appendTextFormatted("%20s of %6d missing %6d", "exclusions", rest_global, -ndiff);
 318         }
 319     }
 320
 321     printMissingInteractionsAtoms(mdlog, cr, top_global, top_local->idef);
 322     write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr, -1, as_rvec_array(x.data()), box);
 323
 324     std::string errorMessage;
 325
 326     if (ndiff_tot > 0)
 327     {
 328         errorMessage =
 329                 "One or more interactions were assigned to multiple domains of the domain "
 330                 "decompostion. Please report this bug.";
 331     }
 332     else
 333     {
 334         errorMessage = formatString(
 335                 "%d of the %d bonded interactions could not be calculated because some atoms "
 336                 "involved moved further apart than the multi-body cut-off distance (%g nm) or the "
 337                 "two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds "
 338                 "also see option -ddcheck",
 339                 -ndiff_tot,
 340                 expectedNumGlobalBondedInteractions,
 341                 dd_cutoff_multibody(dd),
 342                 dd_cutoff_twobody(dd));
 343     }
 344     gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), "%s", errorMessage.c_str());
 345 }
 346
 347 /*! \brief Data to help check local topology construction
 348  *
 349  * Partitioning could incorrectly miss a bonded interaction.
 350  * However, checking for that requires a global communication
 351  * stage, which does not otherwise happen during partitioning. So,
 352  * for performance, we do that alongside the first global energy
 353  * reduction after a new DD is made. These variables handle
 354  * whether the check happens, its input for this domain, output
 355  * across all domains, and the expected value it should match. */
 356 class LocalTopologyChecker::Impl
 357 {
 358 public:
 359     //! Constructor
 360     Impl(const MDLogger& mdlog, const t_commrec* cr, const gmx_mtop_t& mtop, bool useUpdateGroups);
 361
 362     //! Objects used when reporting that interactions are missing
 363     //! {
 364     //! Logger
 365     const MDLogger& mdlog_;
 366     //! Communication record
 367     const t_commrec* cr_;
 368     //! Global system topology
 369     const gmx_mtop_t& mtop_;
 370     //! }
 371
 372     /*! \brief Number of bonded interactions found in the local
 373      * topology for this domain. */
 374     int numBondedInteractionsToReduce_ = 0;
 375     /*! \brief Whether to check at the next global communication
 376      * stage the total number of bonded interactions found.
 377      *
 378      * Cleared after that number is found. */
 379     bool shouldCheckNumberOfBondedInteractions_ = false;
 380     /*! \brief The total number of bonded interactions found in
 381      * the local topology across all domains.
 382      *
 383      * Only has a value after reduction across all ranks, which is
 384      * removed when it is again time to check after a new
 385      * partition. */
 386     std::optional<int> numBondedInteractionsOverAllDomains_;
 387     //! The number of bonded interactions computed from the full system topology
 388     int expectedNumGlobalBondedInteractions_ = 0;
 389 };
 390
 391
 392 /*! \brief Compute the total bonded interaction count
 393  *
 394  * \param[in] mtop              The global system topology
 395  * \param[in] useUpdateGroups   Whether update groups are in use
 396  *
 397  * When using domain decomposition without update groups,
 398  * constraint-type interactions can be split across domains, and so we
 399  * do not consider them in this correctness check. Otherwise, we
 400  * include them.
 401  */
 402 static int computeExpectedNumGlobalBondedInteractions(const gmx_mtop_t& mtop, const bool useUpdateGroups)
 403 {
 404     int expectedNumGlobalBondedInteractions = gmx_mtop_interaction_count(mtop, IF_BOND);
 405     if (useUpdateGroups)
 406     {
 407         expectedNumGlobalBondedInteractions += gmx_mtop_interaction_count(mtop, IF_CONSTRAINT);
 408     }
 409     return expectedNumGlobalBondedInteractions;
 410 }
 411
 412 LocalTopologyChecker::Impl::Impl(const MDLogger&   mdlog,
 413                                  const t_commrec*  cr,
 414                                  const gmx_mtop_t& mtop,
 415                                  const bool        useUpdateGroups) :
 416     mdlog_(mdlog),
 417     cr_(cr),
 418     mtop_(mtop),
 419     expectedNumGlobalBondedInteractions_(computeExpectedNumGlobalBondedInteractions(mtop, useUpdateGroups))
 420 {
 421 }
 422
 423 LocalTopologyChecker::LocalTopologyChecker(const MDLogger&   mdlog,
 424                                            const t_commrec*  cr,
 425                                            const gmx_mtop_t& mtop,
 426                                            const bool        useUpdateGroups) :
 427     impl_(std::make_unique<Impl>(mdlog, cr, mtop, useUpdateGroups))
 428 {
 429 }
 430
 431 LocalTopologyChecker::~LocalTopologyChecker() = default;
 432
 433 LocalTopologyChecker::LocalTopologyChecker(LocalTopologyChecker&&) noexcept = default;
 434
 435 LocalTopologyChecker& LocalTopologyChecker::operator=(LocalTopologyChecker&& other) noexcept
 436 {
 437     impl_ = std::move(other.impl_);
 438     return *this;
 439 }
 440
 441 void LocalTopologyChecker::scheduleCheckOfLocalTopology(const int numBondedInteractionsToReduce)
 442 {
 443     impl_->numBondedInteractionsToReduce_ = numBondedInteractionsToReduce;
 444     // Note that it's possible for this to still be true from the last
 445     // time it was set, e.g. if repartitioning was triggered before
 446     // global communication that would have acted on the true
 447     // value. This could happen for example when replica exchange took
 448     // place soon after a partition.
 449     impl_->shouldCheckNumberOfBondedInteractions_ = true;
 450     // Clear the old global value, which is now invalid
 451     impl_->numBondedInteractionsOverAllDomains_.reset();
 452 }
 453
 454 bool LocalTopologyChecker::shouldCheckNumberOfBondedInteractions() const
 455 {
 456     return impl_->shouldCheckNumberOfBondedInteractions_;
 457 }
 458
 459 int LocalTopologyChecker::numBondedInteractions() const
 460 {
 461     return impl_->numBondedInteractionsToReduce_;
 462 }
 463
 464 void LocalTopologyChecker::setNumberOfBondedInteractionsOverAllDomains(const int newValue)
 465 {
 466     GMX_RELEASE_ASSERT(!impl_->numBondedInteractionsOverAllDomains_.has_value(),
 467                        "Cannot set number of bonded interactions because it is already set");
 468     impl_->numBondedInteractionsOverAllDomains_.emplace(newValue);
 469 }
 470
 471 void LocalTopologyChecker::checkNumberOfBondedInteractions(const gmx_localtop_t* top_local,
 472                                                            ArrayRef<const RVec>  x,
 473                                                            const matrix          box)
 474 {
 475     if (impl_->shouldCheckNumberOfBondedInteractions_)
 476     {
 477         GMX_RELEASE_ASSERT(impl_->numBondedInteractionsOverAllDomains_.has_value(),
 478                            "The check for the total number of bonded interactions requires the "
 479                            "value to have been reduced across all domains");
 480         if (impl_->numBondedInteractionsOverAllDomains_.value() != impl_->expectedNumGlobalBondedInteractions_)
 481         {
 482             dd_print_missing_interactions(impl_->mdlog_,
 483                                           impl_->cr_,
 484                                           impl_->numBondedInteractionsOverAllDomains_.value(),
 485                                           impl_->expectedNumGlobalBondedInteractions_,
 486                                           impl_->mtop_,
 487                                           top_local,
 488                                           x,
 489                                           box); // Does not return
 490         }
 491         // Now that the value is set and the check complete, future
 492         // global communication should not compute the value until
 493         // after the next partitioning.
 494         impl_->shouldCheckNumberOfBondedInteractions_ = false;
 495     }
 496 }
 497
 498 } // namespace gmx