src/gromacs/modularsimulator/computeglobalselement.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,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 /*! \internal \file
  36  * \brief Defines the global reduction element for the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "computeglobalselement.h"
  45
  46 #include "gromacs/domdec/domdec.h"
  47 #include "gromacs/domdec/localtopologychecker.h"
  48 #include "gromacs/gmxlib/network.h"
  49 #include "gromacs/gmxlib/nrnb.h"
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdlib/constr.h"
  52 #include "gromacs/mdlib/md_support.h"
  53 #include "gromacs/mdlib/mdatoms.h"
  54 #include "gromacs/mdlib/stat.h"
  55 #include "gromacs/mdlib/update.h"
  56 #include "gromacs/mdlib/vcm.h"
  57 #include "gromacs/mdtypes/commrec.h"
  58 #include "gromacs/mdtypes/group.h"
  59 #include "gromacs/mdtypes/inputrec.h"
  60 #include "gromacs/mdtypes/mdatom.h"
  61 #include "gromacs/topology/topology.h"
  62
  63 #include "freeenergyperturbationdata.h"
  64 #include "modularsimulator.h"
  65 #include "simulatoralgorithm.h"
  66
  67 namespace gmx
  68 {
  69 template<ComputeGlobalsAlgorithm algorithm>
  70 ComputeGlobalsElement<algorithm>::ComputeGlobalsElement(StatePropagatorData* statePropagatorData,
  71                                                         EnergyData*          energyData,
  72                                                         FreeEnergyPerturbationData* freeEnergyPerturbationData,
  73                                                         SimulationSignals*  signals,
  74                                                         int                 nstglobalcomm,
  75                                                         FILE*               fplog,
  76                                                         const MDLogger&     mdlog,
  77                                                         t_commrec*          cr,
  78                                                         const t_inputrec*   inputrec,
  79                                                         const MDAtoms*      mdAtoms,
  80                                                         t_nrnb*             nrnb,
  81                                                         gmx_wallcycle*      wcycle,
  82                                                         t_forcerec*         fr,
  83                                                         const gmx_mtop_t&   global_top,
  84                                                         Constraints*        constr,
  85                                                         ObservablesReducer* observablesReducer) :
  86     energyReductionStep_(-1),
  87     virialReductionStep_(-1),
  88     vvSchedulingStep_(-1),
  89     doStopCM_(inputrec->comm_mode != ComRemovalAlgorithm::No),
  90     nstcomm_(inputrec->nstcomm),
  91     nstglobalcomm_(nstglobalcomm),
  92     lastStep_(inputrec->nsteps + inputrec->init_step),
  93     initStep_(inputrec->init_step),
  94     nullSignaller_(std::make_unique<SimulationSignaller>(nullptr, nullptr, nullptr, false, false)),
  95     statePropagatorData_(statePropagatorData),
  96     energyData_(energyData),
  97     localTopology_(nullptr),
  98     freeEnergyPerturbationData_(freeEnergyPerturbationData),
  99     vcm_(global_top.groups, *inputrec),
 100     signals_(signals),
 101     fplog_(fplog),
 102     mdlog_(mdlog),
 103     cr_(cr),
 104     inputrec_(inputrec),
 105     top_global_(global_top),
 106     mdAtoms_(mdAtoms),
 107     constr_(constr),
 108     nrnb_(nrnb),
 109     wcycle_(wcycle),
 110     fr_(fr),
 111     observablesReducer_(observablesReducer)
 112 {
 113     reportComRemovalInfo(fplog, vcm_);
 114     gstat_ = global_stat_init(inputrec_);
 115 }
 116
 117 template<ComputeGlobalsAlgorithm algorithm>
 118 ComputeGlobalsElement<algorithm>::~ComputeGlobalsElement()
 119 {
 120     global_stat_destroy(gstat_);
 121 }
 122
 123 template<ComputeGlobalsAlgorithm algorithm>
 124 void ComputeGlobalsElement<algorithm>::elementSetup()
 125 {
 126     GMX_ASSERT(localTopology_, "Setup called before local topology was set.");
 127
 128     if (doStopCM_ && !inputrec_->bContinuation)
 129     {
 130         // To minimize communication, compute_globals computes the COM velocity
 131         // and the kinetic energy for the velocities without COM motion removed.
 132         // Thus to get the kinetic energy without the COM contribution, we need
 133         // to call compute_globals twice.
 134
 135         compute(-1, CGLO_GSTAT | CGLO_STOPCM, nullSignaller_.get(), false, true);
 136
 137         auto v = statePropagatorData_->velocitiesView();
 138         // At initialization, do not pass x with acceleration-correction mode
 139         // to avoid (incorrect) correction of the initial coordinates.
 140         auto x = vcm_.mode == ComRemovalAlgorithm::LinearAccelerationCorrection
 141                          ? ArrayRefWithPadding<RVec>()
 142                          : statePropagatorData_->positionsView();
 143         process_and_stopcm_grp(
 144                 fplog_, &vcm_, *mdAtoms_->mdatoms(), x.unpaddedArrayRef(), v.unpaddedArrayRef());
 145         inc_nrnb(nrnb_, eNR_STOPCM, mdAtoms_->mdatoms()->homenr);
 146     }
 147
 148     unsigned int cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
 149                                | (energyData_->hasReadEkinFromCheckpoint() ? CGLO_READEKIN : 0));
 150
 151     if (algorithm == ComputeGlobalsAlgorithm::VelocityVerlet)
 152     {
 153         cglo_flags |= CGLO_PRESSURE | CGLO_CONSTRAINT;
 154     }
 155
 156     compute(-1, cglo_flags, nullSignaller_.get(), false, true);
 157
 158     // Calculate the initial half step temperature, and save the ekinh_old
 159     for (int i = 0; (i < inputrec_->opts.ngtc); i++)
 160     {
 161         copy_mat(energyData_->ekindata()->tcstat[i].ekinh, energyData_->ekindata()->tcstat[i].ekinh_old);
 162     }
 163 }
 164
 165 template<ComputeGlobalsAlgorithm algorithm>
 166 void ComputeGlobalsElement<algorithm>::scheduleTask(Step                       step,
 167                                                     Time gmx_unused            time,
 168                                                     const RegisterRunFunction& registerRunFunction)
 169 {
 170     const bool needComReduction    = doStopCM_ && do_per_step(step, nstcomm_);
 171     const bool needGlobalReduction = step == energyReductionStep_ || step == virialReductionStep_
 172                                      || needComReduction || do_per_step(step, nstglobalcomm_)
 173                                      || (EI_VV(inputrec_->eI) && inputrecNvtTrotter(inputrec_)
 174                                          && do_per_step(step - 1, nstglobalcomm_));
 175
 176     // TODO: CGLO_GSTAT is only used for needToSumEkinhOld_, i.e. to signal that we do or do not
 177     //       sum the previous kinetic energy. We should simplify / clarify this.
 178
 179     if (algorithm == ComputeGlobalsAlgorithm::LeapFrog)
 180     {
 181         // With Leap-Frog we can skip compute_globals at
 182         // non-communication steps, but we need to calculate
 183         // the kinetic energy one step before communication.
 184
 185         // With leap-frog we also need to compute the half-step
 186         // kinetic energy at the step before we need to write
 187         // the full-step kinetic energy
 188         const bool needEkinAtNextStep = (do_per_step(step + 1, nstglobalcomm_) || step + 1 == lastStep_);
 189
 190         if (!needGlobalReduction && !needEkinAtNextStep)
 191         {
 192             return;
 193         }
 194
 195         const bool doEnergy = step == energyReductionStep_;
 196         int        flags    = (needGlobalReduction ? CGLO_GSTAT : 0) | (doEnergy ? CGLO_ENERGY : 0)
 197                     | (needComReduction ? CGLO_STOPCM : 0) | CGLO_TEMPERATURE | CGLO_PRESSURE
 198                     | CGLO_CONSTRAINT;
 199
 200         // Since we're already communicating at this step, we
 201         // can propagate intra-simulation signals. Note that
 202         // check_nstglobalcomm has the responsibility for
 203         // choosing the value of nstglobalcomm which satisfies
 204         // the need of the different signallers.
 205         const bool doIntraSimSignal = true;
 206         // Disable functionality
 207         const bool doInterSimSignal = false;
 208
 209         // Make signaller to signal stop / reset / checkpointing signals
 210         auto signaller = std::make_shared<SimulationSignaller>(
 211                 signals_, cr_, nullptr, doInterSimSignal, doIntraSimSignal);
 212
 213         registerRunFunction([this, step, flags, signaller = std::move(signaller)]() {
 214             compute(step, flags, signaller.get(), true);
 215         });
 216     }
 217     else if (algorithm == ComputeGlobalsAlgorithm::VelocityVerlet)
 218     {
 219         // For VV, we schedule two calls to compute globals per step.
 220         if (step != vvSchedulingStep_)
 221         {
 222             // This is the first scheduling call for this step (positions & velocities at full time
 223             // step) Set this as the current scheduling step
 224             vvSchedulingStep_ = step;
 225
 226             // For vv, the state at the beginning of the step is positions at time t, velocities at time t - dt/2
 227             // The first velocity propagation (+dt/2) therefore actually corresponds to the previous step.
 228             // So we need information from the last step in the first half of the integration
 229             if (!needGlobalReduction && !do_per_step(step - 1, nstglobalcomm_))
 230             {
 231                 return;
 232             }
 233
 234             const bool doTemperature = step != initStep_ || inputrec_->bContinuation;
 235             const bool doEnergy      = step == energyReductionStep_;
 236
 237             int flags = (needGlobalReduction ? CGLO_GSTAT : 0) | (doEnergy ? CGLO_ENERGY : 0)
 238                         | (doTemperature ? CGLO_TEMPERATURE : 0) | CGLO_PRESSURE | CGLO_CONSTRAINT
 239                         | (needComReduction ? CGLO_STOPCM : 0) | CGLO_SCALEEKIN;
 240
 241             registerRunFunction(
 242                     [this, step, flags]() { compute(step, flags, nullSignaller_.get(), false); });
 243         }
 244         else
 245         {
 246             // second call to compute_globals for this step
 247             // Reset the scheduling step to avoid confusion if scheduling needs
 248             // to be repeated (in case of unexpected simulation termination)
 249             vvSchedulingStep_ = -1;
 250
 251             if (!needGlobalReduction)
 252             {
 253                 return;
 254             }
 255             int flags = CGLO_GSTAT | CGLO_CONSTRAINT;
 256
 257             // Since we're already communicating at this step, we
 258             // can propagate intra-simulation signals. Note that
 259             // check_nstglobalcomm has the responsibility for
 260             // choosing the value of nstglobalcomm which satisfies
 261             // the need of the different signallers.
 262             const bool doIntraSimSignal = true;
 263             // Disable functionality
 264             const bool doInterSimSignal = false;
 265
 266             auto signaller = std::make_shared<SimulationSignaller>(
 267                     signals_, cr_, nullptr, doInterSimSignal, doIntraSimSignal);
 268
 269             registerRunFunction([this, step, flags, signaller = std::move(signaller)]() {
 270                 compute(step, flags, signaller.get(), true);
 271             });
 272         }
 273     }
 274 }
 275
 276 template<ComputeGlobalsAlgorithm algorithm>
 277 void ComputeGlobalsElement<algorithm>::compute(gmx::Step            step,
 278                                                unsigned int         flags,
 279                                                SimulationSignaller* signaller,
 280                                                bool                 useLastBox,
 281                                                bool                 isInit)
 282 {
 283     auto        x       = statePropagatorData_->positionsView().unpaddedArrayRef();
 284     auto        v       = statePropagatorData_->velocitiesView().unpaddedArrayRef();
 285     const auto* box     = statePropagatorData_->constBox();
 286     const auto* lastbox = useLastBox ? statePropagatorData_->constPreviousBox()
 287                                      : statePropagatorData_->constBox();
 288
 289     if (DOMAINDECOMP(cr_) && dd_localTopologyChecker(*cr_->dd).shouldCheckNumberOfBondedInteractions())
 290     {
 291         flags |= CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS;
 292     }
 293     compute_globals(gstat_,
 294                     cr_,
 295                     inputrec_,
 296                     fr_,
 297                     energyData_->ekindata(),
 298                     x,
 299                     v,
 300                     box,
 301                     mdAtoms_->mdatoms(),
 302                     nrnb_,
 303                     &vcm_,
 304                     step != -1 ? wcycle_ : nullptr,
 305                     energyData_->enerdata(),
 306                     energyData_->forceVirial(step),
 307                     energyData_->constraintVirial(step),
 308                     energyData_->totalVirial(step),
 309                     energyData_->pressure(step),
 310                     (((flags & CGLO_ENERGY) != 0) && constr_ != nullptr) ? constr_->rmsdData()
 311                                                                          : gmx::ArrayRef<real>{},
 312                     signaller,
 313                     lastbox,
 314                     energyData_->needToSumEkinhOld(),
 315                     flags,
 316                     step,
 317                     observablesReducer_);
 318     if (DOMAINDECOMP(cr_))
 319     {
 320         dd_localTopologyChecker(cr_->dd)->checkNumberOfBondedInteractions(localTopology_, x, box);
 321     }
 322     if (flags & CGLO_STOPCM && !isInit)
 323     {
 324         process_and_stopcm_grp(fplog_, &vcm_, *mdAtoms_->mdatoms(), x, v);
 325         inc_nrnb(nrnb_, eNR_STOPCM, mdAtoms_->mdatoms()->homenr);
 326     }
 327 }
 328
 329 template<ComputeGlobalsAlgorithm algorithm>
 330 void ComputeGlobalsElement<algorithm>::setTopology(const gmx_localtop_t* top)
 331 {
 332     localTopology_ = top;
 333 }
 334
 335 template<ComputeGlobalsAlgorithm algorithm>
 336 std::optional<SignallerCallback> ComputeGlobalsElement<algorithm>::registerEnergyCallback(EnergySignallerEvent event)
 337 {
 338     if (event == EnergySignallerEvent::EnergyCalculationStep)
 339     {
 340         return [this](Step step, Time /*unused*/) { energyReductionStep_ = step; };
 341     }
 342     if (event == EnergySignallerEvent::VirialCalculationStep)
 343     {
 344         return [this](Step step, Time /*unused*/) { virialReductionStep_ = step; };
 345     }
 346     return std::nullopt;
 347 }
 348
 349 template<ComputeGlobalsAlgorithm algorithm>
 350 std::optional<SignallerCallback>
 351 ComputeGlobalsElement<algorithm>::registerTrajectorySignallerCallback(TrajectoryEvent event)
 352 {
 353     if (event == TrajectoryEvent::EnergyWritingStep)
 354     {
 355         return [this](Step step, Time /*unused*/) { energyReductionStep_ = step; };
 356     }
 357     return std::nullopt;
 358 }
 359
 360 //! Explicit template instantiation
 361 //! \{
 362 template class ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>;
 363 template class ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>;
 364 //! \}
 365
 366 template<>
 367 ISimulatorElement* ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>::getElementPointerImpl(
 368         LegacySimulatorData*                    legacySimulatorData,
 369         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 370         StatePropagatorData*                    statePropagatorData,
 371         EnergyData*                             energyData,
 372         FreeEnergyPerturbationData*             freeEnergyPerturbationData,
 373         GlobalCommunicationHelper*              globalCommunicationHelper,
 374         ObservablesReducer*                     observablesReducer)
 375 {
 376     auto* element = builderHelper->storeElement(
 377             std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>(
 378                     statePropagatorData,
 379                     energyData,
 380                     freeEnergyPerturbationData,
 381                     globalCommunicationHelper->simulationSignals(),
 382                     globalCommunicationHelper->nstglobalcomm(),
 383                     legacySimulatorData->fplog,
 384                     legacySimulatorData->mdlog,
 385                     legacySimulatorData->cr,
 386                     legacySimulatorData->inputrec,
 387                     legacySimulatorData->mdAtoms,
 388                     legacySimulatorData->nrnb,
 389                     legacySimulatorData->wcycle,
 390                     legacySimulatorData->fr,
 391                     legacySimulatorData->top_global,
 392                     legacySimulatorData->constr,
 393                     observablesReducer));
 394
 395     return element;
 396 }
 397
 398 template<>
 399 ISimulatorElement* ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>::getElementPointerImpl(
 400         LegacySimulatorData*                    simulator,
 401         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 402         StatePropagatorData*                    statePropagatorData,
 403         EnergyData*                             energyData,
 404         FreeEnergyPerturbationData*             freeEnergyPerturbationData,
 405         GlobalCommunicationHelper*              globalCommunicationHelper,
 406         ObservablesReducer*                     observablesReducer)
 407 {
 408     // We allow this element to be added multiple times to the call list, but we only want one
 409     // actual element built
 410     static const std::string key("vvComputeGlobalsElement");
 411
 412     const std::optional<std::any> cachedValue = builderHelper->builderData(key);
 413
 414     if (cachedValue)
 415     {
 416         return std::any_cast<ISimulatorElement*>(cachedValue.value());
 417     }
 418     else
 419     {
 420         ISimulatorElement* vvComputeGlobalsElement = builderHelper->storeElement(
 421                 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>(
 422                         statePropagatorData,
 423                         energyData,
 424                         freeEnergyPerturbationData,
 425                         globalCommunicationHelper->simulationSignals(),
 426                         globalCommunicationHelper->nstglobalcomm(),
 427                         simulator->fplog,
 428                         simulator->mdlog,
 429                         simulator->cr,
 430                         simulator->inputrec,
 431                         simulator->mdAtoms,
 432                         simulator->nrnb,
 433                         simulator->wcycle,
 434                         simulator->fr,
 435                         simulator->top_global,
 436                         simulator->constr,
 437                         observablesReducer));
 438         builderHelper->storeBuilderData(key, vvComputeGlobalsElement);
 439         return vvComputeGlobalsElement;
 440     }
 441 }
 442 } // namespace gmx