src/gromacs/modularsimulator/modularsimulator.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020, 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 modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "modularsimulator.h"
  45
  46 #include "gromacs/commandline/filenm.h"
  47 #include "gromacs/domdec/domdec.h"
  48 #include "gromacs/ewald/pme.h"
  49 #include "gromacs/ewald/pme_load_balancing.h"
  50 #include "gromacs/ewald/pme_pp.h"
  51 #include "gromacs/gmxlib/network.h"
  52 #include "gromacs/gmxlib/nrnb.h"
  53 #include "gromacs/listed_forces/listed_forces.h"
  54 #include "gromacs/mdlib/checkpointhandler.h"
  55 #include "gromacs/mdlib/constr.h"
  56 #include "gromacs/mdlib/coupling.h"
  57 #include "gromacs/mdlib/energyoutput.h"
  58 #include "gromacs/mdlib/mdatoms.h"
  59 #include "gromacs/mdlib/resethandler.h"
  60 #include "gromacs/mdlib/update.h"
  61 #include "gromacs/mdrun/replicaexchange.h"
  62 #include "gromacs/mdrun/shellfc.h"
  63 #include "gromacs/mdrunutility/handlerestart.h"
  64 #include "gromacs/mdrunutility/printtime.h"
  65 #include "gromacs/mdtypes/commrec.h"
  66 #include "gromacs/mdtypes/fcdata.h"
  67 #include "gromacs/mdtypes/forcerec.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/mdatom.h"
  70 #include "gromacs/mdtypes/mdrunoptions.h"
  71 #include "gromacs/mdtypes/observableshistory.h"
  72 #include "gromacs/nbnxm/nbnxm.h"
  73 #include "gromacs/topology/mtop_util.h"
  74 #include "gromacs/topology/topology.h"
  75 #include "gromacs/utility/fatalerror.h"
  76
  77 #include "compositesimulatorelement.h"
  78 #include "computeglobalselement.h"
  79 #include "constraintelement.h"
  80 #include "energydata.h"
  81 #include "forceelement.h"
  82 #include "freeenergyperturbationdata.h"
  83 #include "parrinellorahmanbarostat.h"
  84 #include "propagator.h"
  85 #include "signallers.h"
  86 #include "simulatoralgorithm.h"
  87 #include "statepropagatordata.h"
  88 #include "trajectoryelement.h"
  89 #include "vrescalethermostat.h"
  90
  91 namespace gmx
  92 {
  93 void ModularSimulator::run()
  94 {
  95     GMX_LOG(legacySimulatorData_->mdlog.info)
  96             .asParagraph()
  97             .appendText("Using the modular simulator.");
  98
  99     ModularSimulatorAlgorithmBuilder algorithmBuilder(compat::make_not_null(legacySimulatorData_.get()));
 100     auto                             algorithm = algorithmBuilder.build();
 101
 102     while (const auto* task = algorithm.getNextTask())
 103     {
 104         // execute task
 105         (*task)();
 106     }
 107 }
 108
 109 std::unique_ptr<ISimulatorElement> ModularSimulatorAlgorithmBuilder::buildForces(
 110         SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
 111         SignallerBuilder<EnergySignaller>*         energySignallerBuilder,
 112         StatePropagatorData*                       statePropagatorDataPtr,
 113         EnergyData*                                energyDataPtr,
 114         FreeEnergyPerturbationData*                freeEnergyPerturbationDataPtr,
 115         TopologyHolder::Builder*                   topologyHolderBuilder)
 116 {
 117     const bool isVerbose    = legacySimulatorData_->mdrunOptions.verbose;
 118     const bool isDynamicBox = inputrecDynamicBox(legacySimulatorData_->inputrec);
 119
 120     auto forceElement = std::make_unique<ForceElement>(
 121             statePropagatorDataPtr, energyDataPtr, freeEnergyPerturbationDataPtr, isVerbose, isDynamicBox,
 122             legacySimulatorData_->fplog, legacySimulatorData_->cr, legacySimulatorData_->inputrec,
 123             legacySimulatorData_->mdAtoms, legacySimulatorData_->nrnb, legacySimulatorData_->fr,
 124             legacySimulatorData_->wcycle, legacySimulatorData_->runScheduleWork,
 125             legacySimulatorData_->vsite, legacySimulatorData_->imdSession,
 126             legacySimulatorData_->pull_work, legacySimulatorData_->constr,
 127             legacySimulatorData_->top_global, legacySimulatorData_->enforcedRotation);
 128     topologyHolderBuilder->registerClient(forceElement.get());
 129     neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
 130     energySignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
 131
 132     // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 133     return std::move(forceElement);
 134 }
 135
 136 std::unique_ptr<ISimulatorElement> ModularSimulatorAlgorithmBuilder::buildIntegrator(
 137         SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
 138         SignallerBuilder<LastStepSignaller>*       lastStepSignallerBuilder,
 139         SignallerBuilder<EnergySignaller>*         energySignallerBuilder,
 140         SignallerBuilder<LoggingSignaller>*        loggingSignallerBuilder,
 141         SignallerBuilder<TrajectorySignaller>*     trajectorySignallerBuilder,
 142         TrajectoryElementBuilder*                  trajectoryElementBuilder,
 143         std::vector<ICheckpointHelperClient*>*     checkpointClients,
 144         compat::not_null<StatePropagatorData*>     statePropagatorDataPtr,
 145         compat::not_null<EnergyData*>              energyDataPtr,
 146         FreeEnergyPerturbationData*                freeEnergyPerturbationDataPtr,
 147         bool                                       hasReadEkinState,
 148         TopologyHolder::Builder*                   topologyHolderBuilder,
 149         GlobalCommunicationHelper*                 globalCommunicationHelper)
 150 {
 151     auto forceElement = buildForces(neighborSearchSignallerBuilder, energySignallerBuilder,
 152                                     statePropagatorDataPtr, energyDataPtr,
 153                                     freeEnergyPerturbationDataPtr, topologyHolderBuilder);
 154
 155     // list of elements owned by the simulator composite object
 156     std::vector<std::unique_ptr<ISimulatorElement>> elementsOwnershipList;
 157     // call list of the simulator composite object
 158     std::vector<compat::not_null<ISimulatorElement*>> elementCallList;
 159
 160     std::function<void()> needToCheckNumberOfBondedInteractions;
 161     if (legacySimulatorData_->inputrec->eI == eiMD)
 162     {
 163         auto computeGlobalsElement =
 164                 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>(
 165                         statePropagatorDataPtr, energyDataPtr, freeEnergyPerturbationDataPtr,
 166                         globalCommunicationHelper->simulationSignals(),
 167                         globalCommunicationHelper->nstglobalcomm(), legacySimulatorData_->fplog,
 168                         legacySimulatorData_->mdlog, legacySimulatorData_->cr,
 169                         legacySimulatorData_->inputrec, legacySimulatorData_->mdAtoms,
 170                         legacySimulatorData_->nrnb, legacySimulatorData_->wcycle,
 171                         legacySimulatorData_->fr, legacySimulatorData_->top_global,
 172                         legacySimulatorData_->constr, hasReadEkinState);
 173         topologyHolderBuilder->registerClient(computeGlobalsElement.get());
 174         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
 175         trajectorySignallerBuilder->registerSignallerClient(
 176                 compat::make_not_null(computeGlobalsElement.get()));
 177
 178         globalCommunicationHelper->setCheckBondedInteractionsCallback(
 179                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback());
 180
 181         auto propagator = std::make_unique<Propagator<IntegrationStep::LeapFrog>>(
 182                 legacySimulatorData_->inputrec->delta_t, statePropagatorDataPtr,
 183                 legacySimulatorData_->mdAtoms, legacySimulatorData_->wcycle);
 184
 185         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
 186         auto stateElement = compat::make_not_null(statePropagatorDataPtr->element());
 187         trajectoryElementBuilder->registerWriterClient(stateElement);
 188         trajectorySignallerBuilder->registerSignallerClient(stateElement);
 189         lastStepSignallerBuilder->registerSignallerClient(stateElement);
 190         checkpointClients->emplace_back(stateElement);
 191         // we have a full microstate at time t here!
 192         addToCallList(stateElement, elementCallList);
 193         if (legacySimulatorData_->inputrec->etc == etcVRESCALE)
 194         {
 195             // TODO: With increased complexity of the propagator, this will need further development,
 196             //       e.g. using propagators templated for velocity propagation policies and a builder
 197             propagator->setNumVelocityScalingVariables(legacySimulatorData_->inputrec->opts.ngtc);
 198             auto thermostat = std::make_unique<VRescaleThermostat>(
 199                     legacySimulatorData_->inputrec->nsttcouple, -1, false,
 200                     legacySimulatorData_->inputrec->ld_seed, legacySimulatorData_->inputrec->opts.ngtc,
 201                     legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nsttcouple,
 202                     legacySimulatorData_->inputrec->opts.ref_t,
 203                     legacySimulatorData_->inputrec->opts.tau_t, legacySimulatorData_->inputrec->opts.nrdf,
 204                     energyDataPtr, propagator->viewOnVelocityScaling(),
 205                     propagator->velocityScalingCallback(), legacySimulatorData_->state_global,
 206                     legacySimulatorData_->cr, legacySimulatorData_->inputrec->bContinuation);
 207             checkpointClients->emplace_back(thermostat.get());
 208             energyDataPtr->setVRescaleThermostat(thermostat.get());
 209             addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
 210         }
 211
 212         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
 213         if (legacySimulatorData_->inputrec->epc == epcPARRINELLORAHMAN)
 214         {
 215             // Building the PR barostat here since it needs access to the propagator
 216             // and we want to be able to move the propagator object
 217             prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
 218                     legacySimulatorData_->inputrec->nstpcouple, -1,
 219                     legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nstpcouple,
 220                     legacySimulatorData_->inputrec->init_step, propagator->viewOnPRScalingMatrix(),
 221                     propagator->prScalingCallback(), statePropagatorDataPtr, energyDataPtr,
 222                     legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
 223                     legacySimulatorData_->mdAtoms, legacySimulatorData_->state_global,
 224                     legacySimulatorData_->cr, legacySimulatorData_->inputrec->bContinuation);
 225             energyDataPtr->setParrinelloRahamnBarostat(prBarostat.get());
 226             checkpointClients->emplace_back(prBarostat.get());
 227         }
 228         addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
 229         if (legacySimulatorData_->constr)
 230         {
 231             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
 232                     legacySimulatorData_->constr, statePropagatorDataPtr, energyDataPtr,
 233                     freeEnergyPerturbationDataPtr, MASTER(legacySimulatorData_->cr),
 234                     legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
 235                     legacySimulatorData_->mdAtoms->mdatoms());
 236             auto constraintElementPtr = compat::make_not_null(constraintElement.get());
 237             energySignallerBuilder->registerSignallerClient(constraintElementPtr);
 238             trajectorySignallerBuilder->registerSignallerClient(constraintElementPtr);
 239             loggingSignallerBuilder->registerSignallerClient(constraintElementPtr);
 240
 241             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 242         }
 243
 244         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
 245         auto energyElement = compat::make_not_null(energyDataPtr->element());
 246         trajectoryElementBuilder->registerWriterClient(energyElement);
 247         trajectorySignallerBuilder->registerSignallerClient(energyElement);
 248         energySignallerBuilder->registerSignallerClient(energyElement);
 249         checkpointClients->emplace_back(energyElement);
 250         // we have the energies at time t here!
 251         addToCallList(energyElement, elementCallList);
 252         if (prBarostat)
 253         {
 254             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
 255         }
 256     }
 257     else if (legacySimulatorData_->inputrec->eI == eiVV)
 258     {
 259         auto computeGlobalsElement =
 260                 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>(
 261                         statePropagatorDataPtr, energyDataPtr, freeEnergyPerturbationDataPtr,
 262                         globalCommunicationHelper->simulationSignals(),
 263                         globalCommunicationHelper->nstglobalcomm(), legacySimulatorData_->fplog,
 264                         legacySimulatorData_->mdlog, legacySimulatorData_->cr,
 265                         legacySimulatorData_->inputrec, legacySimulatorData_->mdAtoms,
 266                         legacySimulatorData_->nrnb, legacySimulatorData_->wcycle,
 267                         legacySimulatorData_->fr, legacySimulatorData_->top_global,
 268                         legacySimulatorData_->constr, hasReadEkinState);
 269         topologyHolderBuilder->registerClient(computeGlobalsElement.get());
 270         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
 271         trajectorySignallerBuilder->registerSignallerClient(
 272                 compat::make_not_null(computeGlobalsElement.get()));
 273
 274         globalCommunicationHelper->setCheckBondedInteractionsCallback(
 275                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback());
 276
 277         auto propagatorVelocities = std::make_unique<Propagator<IntegrationStep::VelocitiesOnly>>(
 278                 legacySimulatorData_->inputrec->delta_t * 0.5, statePropagatorDataPtr,
 279                 legacySimulatorData_->mdAtoms, legacySimulatorData_->wcycle);
 280         auto propagatorVelocitiesAndPositions =
 281                 std::make_unique<Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>>(
 282                         legacySimulatorData_->inputrec->delta_t, statePropagatorDataPtr,
 283                         legacySimulatorData_->mdAtoms, legacySimulatorData_->wcycle);
 284
 285         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
 286
 287         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
 288         if (legacySimulatorData_->inputrec->epc == epcPARRINELLORAHMAN)
 289         {
 290             // Building the PR barostat here since it needs access to the propagator
 291             // and we want to be able to move the propagator object
 292             prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
 293                     legacySimulatorData_->inputrec->nstpcouple, -1,
 294                     legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nstpcouple,
 295                     legacySimulatorData_->inputrec->init_step,
 296                     propagatorVelocities->viewOnPRScalingMatrix(),
 297                     propagatorVelocities->prScalingCallback(), statePropagatorDataPtr,
 298                     energyDataPtr, legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
 299                     legacySimulatorData_->mdAtoms, legacySimulatorData_->state_global,
 300                     legacySimulatorData_->cr, legacySimulatorData_->inputrec->bContinuation);
 301             energyDataPtr->setParrinelloRahamnBarostat(prBarostat.get());
 302             checkpointClients->emplace_back(prBarostat.get());
 303         }
 304         addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
 305         if (legacySimulatorData_->constr)
 306         {
 307             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Velocities>>(
 308                     legacySimulatorData_->constr, statePropagatorDataPtr, energyDataPtr,
 309                     freeEnergyPerturbationDataPtr, MASTER(legacySimulatorData_->cr),
 310                     legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
 311                     legacySimulatorData_->mdAtoms->mdatoms());
 312             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 313             trajectorySignallerBuilder->registerSignallerClient(
 314                     compat::make_not_null(constraintElement.get()));
 315             loggingSignallerBuilder->registerSignallerClient(
 316                     compat::make_not_null(constraintElement.get()));
 317
 318             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 319         }
 320         addToCallList(compat::make_not_null(computeGlobalsElement.get()), elementCallList);
 321         auto stateElement = compat::make_not_null(statePropagatorDataPtr->element());
 322         trajectoryElementBuilder->registerWriterClient(stateElement);
 323         trajectorySignallerBuilder->registerSignallerClient(stateElement);
 324         lastStepSignallerBuilder->registerSignallerClient(stateElement);
 325         checkpointClients->emplace_back(stateElement);
 326         // we have a full microstate at time t here!
 327         addToCallList(stateElement, elementCallList);
 328         if (legacySimulatorData_->inputrec->etc == etcVRESCALE)
 329         {
 330             // TODO: With increased complexity of the propagator, this will need further development,
 331             //       e.g. using propagators templated for velocity propagation policies and a builder
 332             propagatorVelocitiesAndPositions->setNumVelocityScalingVariables(
 333                     legacySimulatorData_->inputrec->opts.ngtc);
 334             auto thermostat = std::make_unique<VRescaleThermostat>(
 335                     legacySimulatorData_->inputrec->nsttcouple, 0, true,
 336                     legacySimulatorData_->inputrec->ld_seed, legacySimulatorData_->inputrec->opts.ngtc,
 337                     legacySimulatorData_->inputrec->delta_t * legacySimulatorData_->inputrec->nsttcouple,
 338                     legacySimulatorData_->inputrec->opts.ref_t,
 339                     legacySimulatorData_->inputrec->opts.tau_t, legacySimulatorData_->inputrec->opts.nrdf,
 340                     energyDataPtr, propagatorVelocitiesAndPositions->viewOnVelocityScaling(),
 341                     propagatorVelocitiesAndPositions->velocityScalingCallback(),
 342                     legacySimulatorData_->state_global, legacySimulatorData_->cr,
 343                     legacySimulatorData_->inputrec->bContinuation);
 344             checkpointClients->emplace_back(thermostat.get());
 345             energyDataPtr->setVRescaleThermostat(thermostat.get());
 346             addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
 347         }
 348         addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList,
 349                              elementsOwnershipList);
 350         if (legacySimulatorData_->constr)
 351         {
 352             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
 353                     legacySimulatorData_->constr, statePropagatorDataPtr, energyDataPtr,
 354                     freeEnergyPerturbationDataPtr, MASTER(legacySimulatorData_->cr),
 355                     legacySimulatorData_->fplog, legacySimulatorData_->inputrec,
 356                     legacySimulatorData_->mdAtoms->mdatoms());
 357             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 358             trajectorySignallerBuilder->registerSignallerClient(
 359                     compat::make_not_null(constraintElement.get()));
 360             loggingSignallerBuilder->registerSignallerClient(
 361                     compat::make_not_null(constraintElement.get()));
 362
 363             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 364         }
 365         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
 366         auto energyElement = compat::make_not_null(energyDataPtr->element());
 367         trajectoryElementBuilder->registerWriterClient(energyElement);
 368         trajectorySignallerBuilder->registerSignallerClient(energyElement);
 369         energySignallerBuilder->registerSignallerClient(energyElement);
 370         checkpointClients->emplace_back(energyElement);
 371         // we have the energies at time t here!
 372         addToCallList(energyElement, elementCallList);
 373         if (prBarostat)
 374         {
 375             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
 376         }
 377     }
 378     else
 379     {
 380         gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
 381     }
 382
 383     auto integrator = std::make_unique<CompositeSimulatorElement>(std::move(elementCallList),
 384                                                                   std::move(elementsOwnershipList));
 385     // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 386     return std::move(integrator);
 387 }
 388
 389 bool ModularSimulator::isInputCompatible(bool                             exitOnFailure,
 390                                          const t_inputrec*                inputrec,
 391                                          bool                             doRerun,
 392                                          const gmx_mtop_t&                globalTopology,
 393                                          const gmx_multisim_t*            ms,
 394                                          const ReplicaExchangeParameters& replExParams,
 395                                          const t_fcdata*                  fcd,
 396                                          bool                             doEssentialDynamics,
 397                                          bool                             doMembed)
 398 {
 399     auto conditionalAssert = [exitOnFailure](bool condition, const char* message) {
 400         if (exitOnFailure)
 401         {
 402             GMX_RELEASE_ASSERT(condition, message);
 403         }
 404         return condition;
 405     };
 406
 407     bool isInputCompatible = true;
 408
 409     // GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
 410     // such as the leap-frog integrator
 411     const auto modularSimulatorExplicitlyTurnedOn = (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
 412     // GMX_USE_MODULAR_SIMULATOR allows to use disable modular simulator for all uses,
 413     // including the velocity-verlet integrator used by default
 414     const auto modularSimulatorExplicitlyTurnedOff = (getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
 415
 416     GMX_RELEASE_ASSERT(
 417             !(modularSimulatorExplicitlyTurnedOn && modularSimulatorExplicitlyTurnedOff),
 418             "Cannot have both GMX_USE_MODULAR_SIMULATOR=ON and GMX_DISABLE_MODULAR_SIMULATOR=ON. "
 419             "Unset one of the two environment variables to explicitly chose which simulator to "
 420             "use, "
 421             "or unset both to recover default behavior.");
 422
 423     GMX_RELEASE_ASSERT(
 424             !(modularSimulatorExplicitlyTurnedOff && inputrec->eI == eiVV
 425               && inputrec->epc == epcPARRINELLORAHMAN),
 426             "Cannot use a Parrinello-Rahman barostat with md-vv and "
 427             "GMX_DISABLE_MODULAR_SIMULATOR=ON, "
 428             "as the Parrinello-Rahman barostat is not implemented in the legacy simulator. Unset "
 429             "GMX_DISABLE_MODULAR_SIMULATOR or use a different pressure control algorithm.");
 430
 431     isInputCompatible =
 432             isInputCompatible
 433             && conditionalAssert(
 434                        inputrec->eI == eiMD || inputrec->eI == eiVV,
 435                        "Only integrators md and md-vv are supported by the modular simulator.");
 436     isInputCompatible = isInputCompatible
 437                         && conditionalAssert(inputrec->eI != eiMD || modularSimulatorExplicitlyTurnedOn,
 438                                              "Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular "
 439                                              "simulator with integrator md.");
 440     isInputCompatible =
 441             isInputCompatible
 442             && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
 443     isInputCompatible =
 444             isInputCompatible
 445             && conditionalAssert(
 446                        inputrec->etc == etcNO || inputrec->etc == etcVRESCALE,
 447                        "Only v-rescale thermostat is supported by the modular simulator.");
 448     isInputCompatible =
 449             isInputCompatible
 450             && conditionalAssert(
 451                        inputrec->epc == epcNO || inputrec->epc == epcPARRINELLORAHMAN,
 452                        "Only Parrinello-Rahman barostat is supported by the modular simulator.");
 453     isInputCompatible =
 454             isInputCompatible
 455             && conditionalAssert(
 456                        !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec)
 457                          || inputrecNvtTrotter(inputrec)),
 458                        "Legacy Trotter decomposition is not supported by the modular simulator.");
 459     isInputCompatible = isInputCompatible
 460                         && conditionalAssert(inputrec->efep == efepNO || inputrec->efep == efepYES
 461                                                      || inputrec->efep == efepSLOWGROWTH,
 462                                              "Expanded ensemble free energy calculation is not "
 463                                              "supported by the modular simulator.");
 464     isInputCompatible = isInputCompatible
 465                         && conditionalAssert(!inputrec->bPull,
 466                                              "Pulling is not supported by the modular simulator.");
 467     isInputCompatible =
 468             isInputCompatible
 469             && conditionalAssert(inputrec->opts.ngacc == 1 && inputrec->opts.acc[0][XX] == 0.0
 470                                          && inputrec->opts.acc[0][YY] == 0.0
 471                                          && inputrec->opts.acc[0][ZZ] == 0.0 && inputrec->cos_accel == 0.0,
 472                                  "Acceleration is not supported by the modular simulator.");
 473     isInputCompatible =
 474             isInputCompatible
 475             && conditionalAssert(inputrec->opts.ngfrz == 1 && inputrec->opts.nFreeze[0][XX] == 0
 476                                          && inputrec->opts.nFreeze[0][YY] == 0
 477                                          && inputrec->opts.nFreeze[0][ZZ] == 0,
 478                                  "Freeze groups are not supported by the modular simulator.");
 479     isInputCompatible =
 480             isInputCompatible
 481             && conditionalAssert(
 482                        inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
 483                                && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
 484                                && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
 485                                && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
 486                                && inputrec->deform[ZZ][ZZ] == 0.0,
 487                        "Deformation is not supported by the modular simulator.");
 488     isInputCompatible =
 489             isInputCompatible
 490             && conditionalAssert(gmx_mtop_interaction_count(globalTopology, IF_VSITE) == 0,
 491                                  "Virtual sites are not supported by the modular simulator.");
 492     isInputCompatible = isInputCompatible
 493                         && conditionalAssert(!inputrec->bDoAwh,
 494                                              "AWH is not supported by the modular simulator.");
 495     isInputCompatible =
 496             isInputCompatible
 497             && conditionalAssert(gmx_mtop_ftype_count(globalTopology, F_DISRES) == 0,
 498                                  "Distance restraints are not supported by the modular simulator.");
 499     isInputCompatible =
 500             isInputCompatible
 501             && conditionalAssert(
 502                        gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
 503                        "Orientation restraints are not supported by the modular simulator.");
 504     isInputCompatible =
 505             isInputCompatible
 506             && conditionalAssert(ms == nullptr,
 507                                  "Multi-sim are not supported by the modular simulator.");
 508     isInputCompatible =
 509             isInputCompatible
 510             && conditionalAssert(replExParams.exchangeInterval == 0,
 511                                  "Replica exchange is not supported by the modular simulator.");
 512
 513     int numEnsembleRestraintSystems;
 514     if (fcd)
 515     {
 516         numEnsembleRestraintSystems = fcd->disres->nsystems;
 517     }
 518     else
 519     {
 520         auto distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
 521         numEnsembleRestraintSystems =
 522                 (ms != nullptr && distantRestraintEnsembleEnvVar != nullptr)
 523                         ? static_cast<int>(strtol(distantRestraintEnsembleEnvVar, nullptr, 10))
 524                         : 0;
 525     }
 526     isInputCompatible =
 527             isInputCompatible
 528             && conditionalAssert(numEnsembleRestraintSystems <= 1,
 529                                  "Ensemble restraints are not supported by the modular simulator.");
 530     isInputCompatible =
 531             isInputCompatible
 532             && conditionalAssert(!doSimulatedAnnealing(inputrec),
 533                                  "Simulated annealing is not supported by the modular simulator.");
 534     isInputCompatible =
 535             isInputCompatible
 536             && conditionalAssert(!inputrec->bSimTemp,
 537                                  "Simulated tempering is not supported by the modular simulator.");
 538     isInputCompatible = isInputCompatible
 539                         && conditionalAssert(!inputrec->bExpanded,
 540                                              "Expanded ensemble simulations are not supported by "
 541                                              "the modular simulator.");
 542     isInputCompatible =
 543             isInputCompatible
 544             && conditionalAssert(!doEssentialDynamics,
 545                                  "Essential dynamics is not supported by the modular simulator.");
 546     isInputCompatible = isInputCompatible
 547                         && conditionalAssert(inputrec->eSwapCoords == eswapNO,
 548                                              "Ion / water position swapping is not supported by "
 549                                              "the modular simulator.");
 550     isInputCompatible =
 551             isInputCompatible
 552             && conditionalAssert(!inputrec->bIMD,
 553                                  "Interactive MD is not supported by the modular simulator.");
 554     isInputCompatible =
 555             isInputCompatible
 556             && conditionalAssert(!doMembed,
 557                                  "Membrane embedding is not supported by the modular simulator.");
 558     // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
 559     isInputCompatible =
 560             isInputCompatible
 561             && conditionalAssert(
 562                        getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
 563                        "Integration on the GPU is not supported by the modular simulator.");
 564     // Modular simulator is centered around NS updates
 565     // TODO: think how to handle nstlist == 0
 566     isInputCompatible = isInputCompatible
 567                         && conditionalAssert(inputrec->nstlist != 0,
 568                                              "Simulations without neighbor list update are not "
 569                                              "supported by the modular simulator.");
 570     isInputCompatible = isInputCompatible
 571                         && conditionalAssert(!GMX_FAHCORE,
 572                                              "GMX_FAHCORE not supported by the modular simulator.");
 573
 574     return isInputCompatible;
 575 }
 576
 577 ModularSimulator::ModularSimulator(std::unique_ptr<LegacySimulatorData> legacySimulatorData) :
 578     legacySimulatorData_(std::move(legacySimulatorData))
 579 {
 580     checkInputForDisabledFunctionality();
 581 }
 582
 583 void ModularSimulator::checkInputForDisabledFunctionality()
 584 {
 585     isInputCompatible(true, legacySimulatorData_->inputrec,
 586                       legacySimulatorData_->mdrunOptions.rerun, *legacySimulatorData_->top_global,
 587                       legacySimulatorData_->ms, legacySimulatorData_->replExParams,
 588                       &legacySimulatorData_->fr->listedForces->fcdata(),
 589                       opt2bSet("-ei", legacySimulatorData_->nfile, legacySimulatorData_->fnm),
 590                       legacySimulatorData_->membed != nullptr);
 591     if (legacySimulatorData_->observablesHistory->edsamHistory)
 592     {
 593         gmx_fatal(FARGS,
 594                   "The checkpoint is from a run with essential dynamics sampling, "
 595                   "but the current run did not specify the -ei option. "
 596                   "Either specify the -ei option to mdrun, or do not use this checkpoint file.");
 597     }
 598 }
 599
 600 } // namespace gmx