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