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