2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
36 * \brief Defines the modular simulator
38 * \author Pascal Merz <pascal.merz@me.com>
39 * \ingroup module_modularsimulator
44 #include "modularsimulator.h"
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/math/vec.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/forcerec.h"
59 #include "gromacs/mdlib/mdatoms.h"
60 #include "gromacs/mdlib/resethandler.h"
61 #include "gromacs/mdlib/stat.h"
62 #include "gromacs/mdlib/update.h"
63 #include "gromacs/mdrun/replicaexchange.h"
64 #include "gromacs/mdrun/shellfc.h"
65 #include "gromacs/mdrunutility/handlerestart.h"
66 #include "gromacs/mdrunutility/printtime.h"
67 #include "gromacs/mdtypes/commrec.h"
68 #include "gromacs/mdtypes/fcdata.h"
69 #include "gromacs/mdtypes/forcerec.h"
70 #include "gromacs/mdtypes/inputrec.h"
71 #include "gromacs/mdtypes/mdatom.h"
72 #include "gromacs/mdtypes/mdrunoptions.h"
73 #include "gromacs/mdtypes/observableshistory.h"
74 #include "gromacs/mdtypes/state.h"
75 #include "gromacs/nbnxm/nbnxm.h"
76 #include "gromacs/timing/walltime_accounting.h"
77 #include "gromacs/topology/mtop_util.h"
78 #include "gromacs/topology/topology.h"
79 #include "gromacs/utility/cstringutil.h"
80 #include "gromacs/utility/fatalerror.h"
82 #include "compositesimulatorelement.h"
83 #include "computeglobalselement.h"
84 #include "constraintelement.h"
85 #include "energyelement.h"
86 #include "forceelement.h"
87 #include "freeenergyperturbationelement.h"
88 #include "parrinellorahmanbarostat.h"
89 #include "propagator.h"
90 #include "signallers.h"
91 #include "statepropagatordata.h"
92 #include "trajectoryelement.h"
93 #include "vrescalethermostat.h"
97 void ModularSimulator::run()
99 GMX_LOG(mdlog
.info
).asParagraph().appendText("Using the modular simulator.");
100 constructElementsAndSignallers();
102 for (auto& signaller
: signallerCallList_
)
104 signaller
->signallerSetup();
108 domDecHelper_
->setup();
111 for (auto& element
: elementsOwnershipList_
)
113 element
->elementSetup();
115 if (pmeLoadBalanceHelper_
)
117 // State must have been initialized so pmeLoadBalanceHelper_ gets a valid box
118 pmeLoadBalanceHelper_
->setup();
121 while (step_
<= signalHelper_
->lastStep_
)
125 while (!taskQueue_
.empty())
127 auto task
= std::move(taskQueue_
.front());
134 for (auto& element
: elementsOwnershipList_
)
136 element
->elementTeardown();
138 if (pmeLoadBalanceHelper_
)
140 pmeLoadBalanceHelper_
->teardown();
145 void ModularSimulator::simulatorSetup()
147 if (!mdrunOptions
.writeConfout
)
149 // This is on by default, and the main known use case for
150 // turning it off is for convenience in benchmarking, which is
151 // something that should not show up in the general user
156 "The -noconfout functionality is deprecated, and "
157 "may be removed in a future version.");
162 char sbuf
[STEPSTRSIZE
], sbuf2
[STEPSTRSIZE
];
163 std::string timeString
;
164 fprintf(stderr
, "starting mdrun '%s'\n", *(top_global
->name
));
165 if (inputrec
->nsteps
>= 0)
167 timeString
= formatString("%8.1f", static_cast<double>(inputrec
->init_step
+ inputrec
->nsteps
)
168 * inputrec
->delta_t
);
172 timeString
= "infinite";
174 if (inputrec
->init_step
> 0)
176 fprintf(stderr
, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
177 gmx_step_str(inputrec
->init_step
+ inputrec
->nsteps
, sbuf
), timeString
.c_str(),
178 gmx_step_str(inputrec
->init_step
, sbuf2
), inputrec
->init_step
* inputrec
->delta_t
);
182 fprintf(stderr
, "%s steps, %s ps.\n", gmx_step_str(inputrec
->nsteps
, sbuf
),
185 fprintf(fplog
, "\n");
188 walltime_accounting_start_time(walltime_accounting
);
189 wallcycle_start(wcycle
, ewcRUN
);
190 print_start(fplog
, cr
, walltime_accounting
, "mdrun");
192 step_
= inputrec
->init_step
;
195 void ModularSimulator::preStep(Step step
, Time gmx_unused time
, bool isNeighborSearchingStep
)
197 if (stopHandler_
->stoppingAfterCurrentStep(isNeighborSearchingStep
) && step
!= signalHelper_
->lastStep_
)
200 * Stop handler wants to stop after the current step, which was
201 * not known when building the current task queue. This happens
202 * e.g. when a stop is signalled by OS. We therefore want to purge
203 * the task queue now, and re-schedule this step as last step.
206 std::queue
<SimulatorRunFunctionPtr
>().swap(taskQueue_
);
212 resetHandler_
->setSignal(walltime_accounting
);
213 // This is a hack to avoid having to rewrite StopHandler to be a NeighborSearchSignaller
214 // and accept the step as input. Eventually, we want to do that, but currently this would
215 // require introducing NeighborSearchSignaller in the legacy do_md or a lot of code
217 stophandlerIsNSStep_
= isNeighborSearchingStep
;
218 stophandlerCurrentStep_
= step
;
219 stopHandler_
->setSignal();
221 wallcycle_start(wcycle
, ewcSTEP
);
224 void ModularSimulator::postStep(Step step
, Time gmx_unused time
)
229 if (do_per_step(step
, inputrec
->nstlog
))
231 if (fflush(fplog
) != 0)
233 gmx_fatal(FARGS
, "Cannot flush logfile - maybe you are out of disk space?");
237 const bool do_verbose
= mdrunOptions
.verbose
238 && (step
% mdrunOptions
.verboseStepPrintInterval
== 0
239 || step
== inputrec
->init_step
|| step
== signalHelper_
->lastStep_
);
240 // Print the remaining wall clock time for the run
241 if (MASTER(cr
) && (do_verbose
|| gmx_got_usr_signal())
242 && !(pmeLoadBalanceHelper_
&& pmeLoadBalanceHelper_
->pmePrinting()))
244 print_time(stderr
, walltime_accounting
, step
, inputrec
, cr
);
247 double cycles
= wallcycle_stop(wcycle
, ewcSTEP
);
248 if (DOMAINDECOMP(cr
) && wcycle
)
250 dd_cycles_add(cr
->dd
, static_cast<float>(cycles
), ddCyclStep
);
253 resetHandler_
->resetCounters(
254 step
, step
- inputrec
->init_step
, mdlog
, fplog
, cr
, fr
->nbv
.get(), nrnb
, fr
->pmedata
,
255 pmeLoadBalanceHelper_
? pmeLoadBalanceHelper_
->loadBalancingObject() : nullptr, wcycle
,
256 walltime_accounting
);
259 void ModularSimulator::simulatorTeardown()
262 // Stop measuring walltime
263 walltime_accounting_end_time(walltime_accounting
);
265 if (!thisRankHasDuty(cr
, DUTY_PME
))
267 /* Tell the PME only node to finish */
268 gmx_pme_send_finish(cr
);
271 walltime_accounting_set_nsteps_done(walltime_accounting
, step_
- inputrec
->init_step
);
274 void ModularSimulator::populateTaskQueue()
276 auto registerRunFunction
= std::make_unique
<RegisterRunFunction
>(
277 [this](SimulatorRunFunctionPtr ptr
) { taskQueue_
.push(std::move(ptr
)); });
279 Time startTime
= inputrec
->init_t
;
280 Time timeStep
= inputrec
->delta_t
;
281 Time time
= startTime
+ step_
* timeStep
;
283 // Run an initial call to the signallers
284 for (auto& signaller
: signallerCallList_
)
286 signaller
->signal(step_
, time
);
289 if (checkpointHelper_
)
291 checkpointHelper_
->run(step_
, time
);
294 if (pmeLoadBalanceHelper_
)
296 pmeLoadBalanceHelper_
->run(step_
, time
);
300 domDecHelper_
->run(step_
, time
);
305 // local variables for lambda capturing
306 const int step
= step_
;
307 const bool isNSStep
= step
== signalHelper_
->nextNSStep_
;
310 (*registerRunFunction
)(std::make_unique
<SimulatorRunFunction
>(
311 [this, step
, time
, isNSStep
]() { preStep(step
, time
, isNSStep
); }));
312 // register elements for step
313 for (auto& element
: elementCallList_
)
315 element
->scheduleTask(step_
, time
, registerRunFunction
);
317 // register post-step
318 (*registerRunFunction
)(
319 std::make_unique
<SimulatorRunFunction
>([this, step
, time
]() { postStep(step
, time
); }));
323 time
= startTime
+ step_
* timeStep
;
324 for (auto& signaller
: signallerCallList_
)
326 signaller
->signal(step_
, time
);
328 } while (step_
!= signalHelper_
->nextNSStep_
&& step_
<= signalHelper_
->lastStep_
);
331 void ModularSimulator::constructElementsAndSignallers()
333 /* When restarting from a checkpoint, it can be appropriate to
334 * initialize ekind from quantities in the checkpoint. Otherwise,
335 * compute_globals must initialize ekind before the simulation
336 * starts/restarts. However, only the master rank knows what was
337 * found in the checkpoint file, so we have to communicate in
338 * order to coordinate the restart.
340 * TODO (modular) This should become obsolete when checkpoint reading
341 * happens within the modular simulator framework: The energy
342 * element should read its data from the checkpoint file pointer,
343 * and signal to the compute globals element if it needs anything
346 * TODO (legacy) Consider removing this communication if/when checkpoint
347 * reading directly follows .tpr reading, because all ranks can
348 * agree on hasReadEkinState at that time.
350 bool hasReadEkinState
= MASTER(cr
) ? state_global
->ekinstate
.hasReadEkinState
: false;
353 gmx_bcast(sizeof(hasReadEkinState
), &hasReadEkinState
, cr
->mpi_comm_mygroup
);
355 if (hasReadEkinState
)
357 restore_ekinstate_from_state(cr
, ekind
, &state_global
->ekinstate
);
361 * Build data structures
364 std::make_unique
<TopologyHolder
>(*top_global
, cr
, inputrec
, fr
, mdAtoms
, constr
, vsite
);
366 std::unique_ptr
<FreeEnergyPerturbationElement
> freeEnergyPerturbationElement
= nullptr;
367 FreeEnergyPerturbationElement
* freeEnergyPerturbationElementPtr
= nullptr;
368 if (inputrec
->efep
!= efepNO
)
370 freeEnergyPerturbationElement
=
371 std::make_unique
<FreeEnergyPerturbationElement
>(fplog
, inputrec
, mdAtoms
);
372 freeEnergyPerturbationElementPtr
= freeEnergyPerturbationElement
.get();
375 auto statePropagatorData
= std::make_unique
<StatePropagatorData
>(
376 top_global
->natoms
, fplog
, cr
, state_global
, inputrec
->nstxout
, inputrec
->nstvout
,
377 inputrec
->nstfout
, inputrec
->nstxout_compressed
, fr
->nbv
->useGpu(),
378 freeEnergyPerturbationElementPtr
, topologyHolder_
.get(), fr
->bMolPBC
,
379 mdrunOptions
.writeConfout
, opt2fn("-c", nfile
, fnm
), inputrec
, mdAtoms
->mdatoms());
380 auto statePropagatorDataPtr
= compat::make_not_null(statePropagatorData
.get());
382 auto energyElement
= std::make_unique
<EnergyElement
>(
383 statePropagatorDataPtr
, freeEnergyPerturbationElementPtr
, top_global
, inputrec
, mdAtoms
,
384 enerd
, ekind
, constr
, fplog
, fcd
, mdModulesNotifier
, MASTER(cr
), observablesHistory
,
386 auto energyElementPtr
= compat::make_not_null(energyElement
.get());
391 const bool simulationsShareState
= false;
392 stopHandler_
= stopHandlerBuilder
->getStopHandlerMD(
393 compat::not_null
<SimulationSignal
*>(&signals_
[eglsSTOPCOND
]), simulationsShareState
,
394 MASTER(cr
), inputrec
->nstlist
, mdrunOptions
.reproducible
, nstglobalcomm_
,
395 mdrunOptions
.maximumHoursToRun
, inputrec
->nstlist
== 0, fplog
, stophandlerCurrentStep_
,
396 stophandlerIsNSStep_
, walltime_accounting
);
399 * Create simulator builders
401 SignallerBuilder
<NeighborSearchSignaller
> neighborSearchSignallerBuilder
;
402 SignallerBuilder
<LastStepSignaller
> lastStepSignallerBuilder
;
403 SignallerBuilder
<LoggingSignaller
> loggingSignallerBuilder
;
404 SignallerBuilder
<EnergySignaller
> energySignallerBuilder
;
405 TrajectoryElementBuilder trajectoryElementBuilder
;
408 * Register data structures to signallers
410 trajectoryElementBuilder
.registerWriterClient(statePropagatorDataPtr
);
411 trajectoryElementBuilder
.registerSignallerClient(statePropagatorDataPtr
);
412 lastStepSignallerBuilder
.registerSignallerClient(statePropagatorDataPtr
);
414 trajectoryElementBuilder
.registerWriterClient(energyElementPtr
);
415 trajectoryElementBuilder
.registerSignallerClient(energyElementPtr
);
416 energySignallerBuilder
.registerSignallerClient(energyElementPtr
);
418 // Register the simulator itself to the neighbor search / last step signaller
419 neighborSearchSignallerBuilder
.registerSignallerClient(compat::make_not_null(signalHelper_
.get()));
420 lastStepSignallerBuilder
.registerSignallerClient(compat::make_not_null(signalHelper_
.get()));
423 * Build integrator - this takes care of force calculation, propagation,
424 * constraining, and of the place the statePropagatorData and the energy element
425 * have a full timestep state.
427 // TODO: Make a CheckpointHelperBuilder
428 std::vector
<ICheckpointHelperClient
*> checkpointClients
= { statePropagatorDataPtr
, energyElementPtr
,
429 freeEnergyPerturbationElementPtr
};
430 CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback
= nullptr;
432 buildIntegrator(&neighborSearchSignallerBuilder
, &energySignallerBuilder
,
433 &loggingSignallerBuilder
, &trajectoryElementBuilder
, &checkpointClients
,
434 &checkBondedInteractionsCallback
, statePropagatorDataPtr
,
435 energyElementPtr
, freeEnergyPerturbationElementPtr
, hasReadEkinState
);
438 * Build infrastructure elements
441 if (PmeLoadBalanceHelper::doPmeLoadBalancing(mdrunOptions
, inputrec
, fr
))
443 pmeLoadBalanceHelper_
= std::make_unique
<PmeLoadBalanceHelper
>(
444 mdrunOptions
.verbose
, statePropagatorDataPtr
, fplog
, cr
, mdlog
, inputrec
, wcycle
, fr
);
445 neighborSearchSignallerBuilder
.registerSignallerClient(
446 compat::make_not_null(pmeLoadBalanceHelper_
.get()));
449 if (DOMAINDECOMP(cr
))
451 GMX_ASSERT(checkBondedInteractionsCallback
,
452 "Domain decomposition needs a callback for check the number of bonded "
454 domDecHelper_
= std::make_unique
<DomDecHelper
>(
455 mdrunOptions
.verbose
, mdrunOptions
.verboseStepPrintInterval
, statePropagatorDataPtr
,
456 topologyHolder_
.get(), std::move(checkBondedInteractionsCallback
), nstglobalcomm_
, fplog
,
457 cr
, mdlog
, constr
, inputrec
, mdAtoms
, nrnb
, wcycle
, fr
, vsite
, imdSession
, pull_work
);
458 neighborSearchSignallerBuilder
.registerSignallerClient(compat::make_not_null(domDecHelper_
.get()));
461 const bool simulationsShareResetCounters
= false;
462 resetHandler_
= std::make_unique
<ResetHandler
>(
463 compat::make_not_null
<SimulationSignal
*>(&signals_
[eglsRESETCOUNTERS
]),
464 simulationsShareResetCounters
, inputrec
->nsteps
, MASTER(cr
),
465 mdrunOptions
.timingOptions
.resetHalfway
, mdrunOptions
.maximumHoursToRun
, mdlog
, wcycle
,
466 walltime_accounting
);
469 * Build signaller list
471 * Note that as signallers depend on each others, the order of calling the signallers
472 * matters. It is the responsibility of this builder to ensure that the order is
475 auto energySignaller
= energySignallerBuilder
.build(
476 inputrec
->nstcalcenergy
, inputrec
->fepvals
->nstdhdl
, inputrec
->nstpcouple
);
477 trajectoryElementBuilder
.registerSignallerClient(compat::make_not_null(energySignaller
.get()));
478 loggingSignallerBuilder
.registerSignallerClient(compat::make_not_null(energySignaller
.get()));
479 auto trajectoryElement
= trajectoryElementBuilder
.build(
480 fplog
, nfile
, fnm
, mdrunOptions
, cr
, outputProvider
, mdModulesNotifier
, inputrec
,
481 top_global
, oenv
, wcycle
, startingBehavior
, simulationsShareState
);
482 loggingSignallerBuilder
.registerSignallerClient(compat::make_not_null(trajectoryElement
.get()));
484 // Add checkpoint helper here since we need a pointer to the trajectory element and
485 // need to register it with the lastStepSignallerBuilder
486 auto checkpointHandler
= std::make_unique
<CheckpointHandler
>(
487 compat::make_not_null
<SimulationSignal
*>(&signals_
[eglsCHKPT
]), simulationsShareState
,
488 inputrec
->nstlist
== 0, MASTER(cr
), mdrunOptions
.writeConfout
,
489 mdrunOptions
.checkpointOptions
.period
);
490 checkpointHelper_
= std::make_unique
<CheckpointHelper
>(
491 std::move(checkpointClients
), std::move(checkpointHandler
), inputrec
->init_step
,
492 trajectoryElement
.get(), top_global
->natoms
, fplog
, cr
, observablesHistory
,
493 walltime_accounting
, state_global
, mdrunOptions
.writeConfout
);
494 lastStepSignallerBuilder
.registerSignallerClient(compat::make_not_null(checkpointHelper_
.get()));
496 lastStepSignallerBuilder
.registerSignallerClient(compat::make_not_null(trajectoryElement
.get()));
497 auto loggingSignaller
=
498 loggingSignallerBuilder
.build(inputrec
->nstlog
, inputrec
->init_step
, inputrec
->init_t
);
499 lastStepSignallerBuilder
.registerSignallerClient(compat::make_not_null(loggingSignaller
.get()));
500 auto lastStepSignaller
=
501 lastStepSignallerBuilder
.build(inputrec
->nsteps
, inputrec
->init_step
, stopHandler_
.get());
502 neighborSearchSignallerBuilder
.registerSignallerClient(compat::make_not_null(lastStepSignaller
.get()));
503 auto neighborSearchSignaller
= neighborSearchSignallerBuilder
.build(
504 inputrec
->nstlist
, inputrec
->init_step
, inputrec
->init_t
);
506 addToCallListAndMove(std::move(neighborSearchSignaller
), signallerCallList_
, signallersOwnershipList_
);
507 addToCallListAndMove(std::move(lastStepSignaller
), signallerCallList_
, signallersOwnershipList_
);
508 addToCallListAndMove(std::move(loggingSignaller
), signallerCallList_
, signallersOwnershipList_
);
509 addToCallList(trajectoryElement
, signallerCallList_
);
510 addToCallListAndMove(std::move(energySignaller
), signallerCallList_
, signallersOwnershipList_
);
513 * Build the element list
515 * This is the actual sequence of (non-infrastructure) elements to be run.
516 * For NVE, the trajectory element is used outside of the integrator
517 * (composite) element, as well as the checkpoint helper. The checkpoint
518 * helper should be on top of the loop, and is only part of the simulator
519 * call list to be able to react to the last step being signalled.
521 addToCallList(checkpointHelper_
, elementCallList_
);
522 if (freeEnergyPerturbationElement
)
524 addToCallListAndMove(std::move(freeEnergyPerturbationElement
), elementCallList_
,
525 elementsOwnershipList_
);
527 addToCallListAndMove(std::move(integrator
), elementCallList_
, elementsOwnershipList_
);
528 addToCallListAndMove(std::move(trajectoryElement
), elementCallList_
, elementsOwnershipList_
);
529 // for vv, we need to setup statePropagatorData after the compute
530 // globals so that we reset the right velocities
531 // TODO: Avoid this by getting rid of the need of resetting velocities in vv
532 elementsOwnershipList_
.emplace_back(std::move(statePropagatorData
));
533 elementsOwnershipList_
.emplace_back(std::move(energyElement
));
536 std::unique_ptr
<ISimulatorElement
>
537 ModularSimulator::buildForces(SignallerBuilder
<NeighborSearchSignaller
>* neighborSearchSignallerBuilder
,
538 SignallerBuilder
<EnergySignaller
>* energySignallerBuilder
,
539 StatePropagatorData
* statePropagatorDataPtr
,
540 EnergyElement
* energyElementPtr
,
541 FreeEnergyPerturbationElement
* freeEnergyPerturbationElement
)
543 const bool isVerbose
= mdrunOptions
.verbose
;
544 const bool isDynamicBox
= inputrecDynamicBox(inputrec
);
546 auto forceElement
= std::make_unique
<ForceElement
>(
547 statePropagatorDataPtr
, energyElementPtr
, freeEnergyPerturbationElement
, isVerbose
,
548 isDynamicBox
, fplog
, cr
, inputrec
, mdAtoms
, nrnb
, fr
, fcd
, wcycle
, runScheduleWork
, vsite
,
549 imdSession
, pull_work
, constr
, &topologyHolder_
->globalTopology(), enforcedRotation
);
550 topologyHolder_
->registerClient(forceElement
.get());
551 neighborSearchSignallerBuilder
->registerSignallerClient(compat::make_not_null(forceElement
.get()));
552 energySignallerBuilder
->registerSignallerClient(compat::make_not_null(forceElement
.get()));
554 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
555 return std::move(forceElement
);
558 std::unique_ptr
<ISimulatorElement
> ModularSimulator::buildIntegrator(
559 SignallerBuilder
<NeighborSearchSignaller
>* neighborSearchSignallerBuilder
,
560 SignallerBuilder
<EnergySignaller
>* energySignallerBuilder
,
561 SignallerBuilder
<LoggingSignaller
>* loggingSignallerBuilder
,
562 TrajectoryElementBuilder
* trajectoryElementBuilder
,
563 std::vector
<ICheckpointHelperClient
*>* checkpointClients
,
564 CheckBondedInteractionsCallbackPtr
* checkBondedInteractionsCallback
,
565 compat::not_null
<StatePropagatorData
*> statePropagatorDataPtr
,
566 compat::not_null
<EnergyElement
*> energyElementPtr
,
567 FreeEnergyPerturbationElement
* freeEnergyPerturbationElementPtr
,
568 bool hasReadEkinState
)
571 buildForces(neighborSearchSignallerBuilder
, energySignallerBuilder
,
572 statePropagatorDataPtr
, energyElementPtr
, freeEnergyPerturbationElementPtr
);
574 // list of elements owned by the simulator composite object
575 std::vector
<std::unique_ptr
<ISimulatorElement
>> elementsOwnershipList
;
576 // call list of the simulator composite object
577 std::vector
<compat::not_null
<ISimulatorElement
*>> elementCallList
;
579 std::function
<void()> needToCheckNumberOfBondedInteractions
;
580 if (inputrec
->eI
== eiMD
)
582 auto computeGlobalsElement
=
583 std::make_unique
<ComputeGlobalsElement
<ComputeGlobalsAlgorithm::LeapFrog
>>(
584 statePropagatorDataPtr
, energyElementPtr
, freeEnergyPerturbationElementPtr
,
585 &signals_
, nstglobalcomm_
, fplog
, mdlog
, cr
, inputrec
, mdAtoms
, nrnb
,
586 wcycle
, fr
, &topologyHolder_
->globalTopology(), constr
, hasReadEkinState
);
587 topologyHolder_
->registerClient(computeGlobalsElement
.get());
588 energySignallerBuilder
->registerSignallerClient(compat::make_not_null(computeGlobalsElement
.get()));
589 trajectoryElementBuilder
->registerSignallerClient(
590 compat::make_not_null(computeGlobalsElement
.get()));
592 *checkBondedInteractionsCallback
=
593 computeGlobalsElement
->getCheckNumberOfBondedInteractionsCallback();
595 auto propagator
= std::make_unique
<Propagator
<IntegrationStep::LeapFrog
>>(
596 inputrec
->delta_t
, statePropagatorDataPtr
, mdAtoms
, wcycle
);
598 addToCallListAndMove(std::move(forceElement
), elementCallList
, elementsOwnershipList
);
599 addToCallList(statePropagatorDataPtr
, elementCallList
); // we have a full microstate at time t here!
600 if (inputrec
->etc
== etcVRESCALE
)
602 // TODO: With increased complexity of the propagator, this will need further development,
603 // e.g. using propagators templated for velocity propagation policies and a builder
604 propagator
->setNumVelocityScalingVariables(inputrec
->opts
.ngtc
);
605 auto thermostat
= std::make_unique
<VRescaleThermostat
>(
606 inputrec
->nsttcouple
, -1, false, inputrec
->ld_seed
, inputrec
->opts
.ngtc
,
607 inputrec
->delta_t
* inputrec
->nsttcouple
, inputrec
->opts
.ref_t
, inputrec
->opts
.tau_t
,
608 inputrec
->opts
.nrdf
, energyElementPtr
, propagator
->viewOnVelocityScaling(),
609 propagator
->velocityScalingCallback(), state_global
, cr
, inputrec
->bContinuation
);
610 checkpointClients
->emplace_back(thermostat
.get());
611 energyElementPtr
->setVRescaleThermostat(thermostat
.get());
612 addToCallListAndMove(std::move(thermostat
), elementCallList
, elementsOwnershipList
);
615 std::unique_ptr
<ParrinelloRahmanBarostat
> prBarostat
= nullptr;
616 if (inputrec
->epc
== epcPARRINELLORAHMAN
)
618 // Building the PR barostat here since it needs access to the propagator
619 // and we want to be able to move the propagator object
620 prBarostat
= std::make_unique
<ParrinelloRahmanBarostat
>(
621 inputrec
->nstpcouple
, -1, inputrec
->delta_t
* inputrec
->nstpcouple
,
622 inputrec
->init_step
, propagator
->viewOnPRScalingMatrix(),
623 propagator
->prScalingCallback(), statePropagatorDataPtr
, energyElementPtr
,
624 fplog
, inputrec
, mdAtoms
, state_global
, cr
, inputrec
->bContinuation
);
625 energyElementPtr
->setParrinelloRahamnBarostat(prBarostat
.get());
626 checkpointClients
->emplace_back(prBarostat
.get());
628 addToCallListAndMove(std::move(propagator
), elementCallList
, elementsOwnershipList
);
631 auto constraintElement
= std::make_unique
<ConstraintsElement
<ConstraintVariable::Positions
>>(
632 constr
, statePropagatorDataPtr
, energyElementPtr
, freeEnergyPerturbationElementPtr
,
633 MASTER(cr
), fplog
, inputrec
, mdAtoms
->mdatoms());
634 auto constraintElementPtr
= compat::make_not_null(constraintElement
.get());
635 energySignallerBuilder
->registerSignallerClient(constraintElementPtr
);
636 trajectoryElementBuilder
->registerSignallerClient(constraintElementPtr
);
637 loggingSignallerBuilder
->registerSignallerClient(constraintElementPtr
);
639 addToCallListAndMove(std::move(constraintElement
), elementCallList
, elementsOwnershipList
);
642 addToCallListAndMove(std::move(computeGlobalsElement
), elementCallList
, elementsOwnershipList
);
643 addToCallList(energyElementPtr
, elementCallList
); // we have the energies at time t here!
646 addToCallListAndMove(std::move(prBarostat
), elementCallList
, elementsOwnershipList
);
649 else if (inputrec
->eI
== eiVV
)
651 auto computeGlobalsElement
=
652 std::make_unique
<ComputeGlobalsElement
<ComputeGlobalsAlgorithm::VelocityVerlet
>>(
653 statePropagatorDataPtr
, energyElementPtr
, freeEnergyPerturbationElementPtr
,
654 &signals_
, nstglobalcomm_
, fplog
, mdlog
, cr
, inputrec
, mdAtoms
, nrnb
,
655 wcycle
, fr
, &topologyHolder_
->globalTopology(), constr
, hasReadEkinState
);
656 topologyHolder_
->registerClient(computeGlobalsElement
.get());
657 energySignallerBuilder
->registerSignallerClient(compat::make_not_null(computeGlobalsElement
.get()));
658 trajectoryElementBuilder
->registerSignallerClient(
659 compat::make_not_null(computeGlobalsElement
.get()));
661 *checkBondedInteractionsCallback
=
662 computeGlobalsElement
->getCheckNumberOfBondedInteractionsCallback();
664 auto propagatorVelocities
= std::make_unique
<Propagator
<IntegrationStep::VelocitiesOnly
>>(
665 inputrec
->delta_t
* 0.5, statePropagatorDataPtr
, mdAtoms
, wcycle
);
666 auto propagatorVelocitiesAndPositions
=
667 std::make_unique
<Propagator
<IntegrationStep::VelocityVerletPositionsAndVelocities
>>(
668 inputrec
->delta_t
, statePropagatorDataPtr
, mdAtoms
, wcycle
);
670 addToCallListAndMove(std::move(forceElement
), elementCallList
, elementsOwnershipList
);
672 std::unique_ptr
<ParrinelloRahmanBarostat
> prBarostat
= nullptr;
673 if (inputrec
->epc
== epcPARRINELLORAHMAN
)
675 // Building the PR barostat here since it needs access to the propagator
676 // and we want to be able to move the propagator object
677 prBarostat
= std::make_unique
<ParrinelloRahmanBarostat
>(
678 inputrec
->nstpcouple
, -1, inputrec
->delta_t
* inputrec
->nstpcouple
,
679 inputrec
->init_step
, propagatorVelocities
->viewOnPRScalingMatrix(),
680 propagatorVelocities
->prScalingCallback(), statePropagatorDataPtr
, energyElementPtr
,
681 fplog
, inputrec
, mdAtoms
, state_global
, cr
, inputrec
->bContinuation
);
682 energyElementPtr
->setParrinelloRahamnBarostat(prBarostat
.get());
683 checkpointClients
->emplace_back(prBarostat
.get());
685 addToCallListAndMove(std::move(propagatorVelocities
), elementCallList
, elementsOwnershipList
);
688 auto constraintElement
= std::make_unique
<ConstraintsElement
<ConstraintVariable::Velocities
>>(
689 constr
, statePropagatorDataPtr
, energyElementPtr
, freeEnergyPerturbationElementPtr
,
690 MASTER(cr
), fplog
, inputrec
, mdAtoms
->mdatoms());
691 energySignallerBuilder
->registerSignallerClient(compat::make_not_null(constraintElement
.get()));
692 trajectoryElementBuilder
->registerSignallerClient(
693 compat::make_not_null(constraintElement
.get()));
694 loggingSignallerBuilder
->registerSignallerClient(
695 compat::make_not_null(constraintElement
.get()));
697 addToCallListAndMove(std::move(constraintElement
), elementCallList
, elementsOwnershipList
);
699 addToCallList(compat::make_not_null(computeGlobalsElement
.get()), elementCallList
);
700 addToCallList(statePropagatorDataPtr
, elementCallList
); // we have a full microstate at time t here!
701 if (inputrec
->etc
== etcVRESCALE
)
703 // TODO: With increased complexity of the propagator, this will need further development,
704 // e.g. using propagators templated for velocity propagation policies and a builder
705 propagatorVelocitiesAndPositions
->setNumVelocityScalingVariables(inputrec
->opts
.ngtc
);
706 auto thermostat
= std::make_unique
<VRescaleThermostat
>(
707 inputrec
->nsttcouple
, 0, true, inputrec
->ld_seed
, inputrec
->opts
.ngtc
,
708 inputrec
->delta_t
* inputrec
->nsttcouple
, inputrec
->opts
.ref_t
,
709 inputrec
->opts
.tau_t
, inputrec
->opts
.nrdf
, energyElementPtr
,
710 propagatorVelocitiesAndPositions
->viewOnVelocityScaling(),
711 propagatorVelocitiesAndPositions
->velocityScalingCallback(), state_global
, cr
,
712 inputrec
->bContinuation
);
713 checkpointClients
->emplace_back(thermostat
.get());
714 energyElementPtr
->setVRescaleThermostat(thermostat
.get());
715 addToCallListAndMove(std::move(thermostat
), elementCallList
, elementsOwnershipList
);
717 addToCallListAndMove(std::move(propagatorVelocitiesAndPositions
), elementCallList
,
718 elementsOwnershipList
);
721 auto constraintElement
= std::make_unique
<ConstraintsElement
<ConstraintVariable::Positions
>>(
722 constr
, statePropagatorDataPtr
, energyElementPtr
, freeEnergyPerturbationElementPtr
,
723 MASTER(cr
), fplog
, inputrec
, mdAtoms
->mdatoms());
724 energySignallerBuilder
->registerSignallerClient(compat::make_not_null(constraintElement
.get()));
725 trajectoryElementBuilder
->registerSignallerClient(
726 compat::make_not_null(constraintElement
.get()));
727 loggingSignallerBuilder
->registerSignallerClient(
728 compat::make_not_null(constraintElement
.get()));
730 addToCallListAndMove(std::move(constraintElement
), elementCallList
, elementsOwnershipList
);
732 addToCallListAndMove(std::move(computeGlobalsElement
), elementCallList
, elementsOwnershipList
);
733 addToCallList(energyElementPtr
, elementCallList
); // we have the energies at time t here!
736 addToCallListAndMove(std::move(prBarostat
), elementCallList
, elementsOwnershipList
);
741 gmx_fatal(FARGS
, "Integrator not implemented for the modular simulator.");
744 auto integrator
= std::make_unique
<CompositeSimulatorElement
>(std::move(elementCallList
),
745 std::move(elementsOwnershipList
));
746 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
747 return std::move(integrator
);
750 bool ModularSimulator::isInputCompatible(bool exitOnFailure
,
751 const t_inputrec
* inputrec
,
753 const gmx_mtop_t
& globalTopology
,
754 const gmx_multisim_t
* ms
,
755 const ReplicaExchangeParameters
& replExParams
,
757 bool doEssentialDynamics
,
760 auto conditionalAssert
= [exitOnFailure
](bool condition
, const char* message
) {
763 GMX_RELEASE_ASSERT(condition
, message
);
768 bool isInputCompatible
= true;
770 // GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
771 // such as the leap-frog integrator
772 const auto modularSimulatorExplicitlyTurnedOn
= (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
773 // GMX_USE_MODULAR_SIMULATOR allows to use disable modular simulator for all uses,
774 // including the velocity-verlet integrator used by default
775 const auto modularSimulatorExplicitlyTurnedOff
= (getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
778 !(modularSimulatorExplicitlyTurnedOn
&& modularSimulatorExplicitlyTurnedOff
),
779 "Cannot have both GMX_USE_MODULAR_SIMULATOR=ON and GMX_DISABLE_MODULAR_SIMULATOR=ON. "
780 "Unset one of the two environment variables to explicitly chose which simulator to "
782 "or unset both to recover default behavior.");
785 !(modularSimulatorExplicitlyTurnedOff
&& inputrec
->eI
== eiVV
786 && inputrec
->epc
== epcPARRINELLORAHMAN
),
787 "Cannot use a Parrinello-Rahman barostat with md-vv and "
788 "GMX_DISABLE_MODULAR_SIMULATOR=ON, "
789 "as the Parrinello-Rahman barostat is not implemented in the legacy simulator. Unset "
790 "GMX_DISABLE_MODULAR_SIMULATOR or use a different pressure control algorithm.");
794 && conditionalAssert(
795 inputrec
->eI
== eiMD
|| inputrec
->eI
== eiVV
,
796 "Only integrators md and md-vv are supported by the modular simulator.");
797 isInputCompatible
= isInputCompatible
798 && conditionalAssert(inputrec
->eI
!= eiMD
|| modularSimulatorExplicitlyTurnedOn
,
799 "Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular "
800 "simulator with integrator md.");
803 && conditionalAssert(!doRerun
, "Rerun is not supported by the modular simulator.");
806 && conditionalAssert(
807 inputrec
->etc
== etcNO
|| inputrec
->etc
== etcVRESCALE
,
808 "Only v-rescale thermostat is supported by the modular simulator.");
811 && conditionalAssert(
812 inputrec
->epc
== epcNO
|| inputrec
->epc
== epcPARRINELLORAHMAN
,
813 "Only Parrinello-Rahman barostat is supported by the modular simulator.");
816 && conditionalAssert(
817 !(inputrecNptTrotter(inputrec
) || inputrecNphTrotter(inputrec
)
818 || inputrecNvtTrotter(inputrec
)),
819 "Legacy Trotter decomposition is not supported by the modular simulator.");
820 isInputCompatible
= isInputCompatible
821 && conditionalAssert(inputrec
->efep
== efepNO
|| inputrec
->efep
== efepYES
822 || inputrec
->efep
== efepSLOWGROWTH
,
823 "Expanded ensemble free energy calculation is not "
824 "supported by the modular simulator.");
825 isInputCompatible
= isInputCompatible
826 && conditionalAssert(!inputrec
->bPull
,
827 "Pulling is not supported by the modular simulator.");
830 && conditionalAssert(inputrec
->opts
.ngacc
== 1 && inputrec
->opts
.acc
[0][XX
] == 0.0
831 && inputrec
->opts
.acc
[0][YY
] == 0.0
832 && inputrec
->opts
.acc
[0][ZZ
] == 0.0 && inputrec
->cos_accel
== 0.0,
833 "Acceleration is not supported by the modular simulator.");
836 && conditionalAssert(inputrec
->opts
.ngfrz
== 1 && inputrec
->opts
.nFreeze
[0][XX
] == 0
837 && inputrec
->opts
.nFreeze
[0][YY
] == 0
838 && inputrec
->opts
.nFreeze
[0][ZZ
] == 0,
839 "Freeze groups are not supported by the modular simulator.");
842 && conditionalAssert(
843 inputrec
->deform
[XX
][XX
] == 0.0 && inputrec
->deform
[XX
][YY
] == 0.0
844 && inputrec
->deform
[XX
][ZZ
] == 0.0 && inputrec
->deform
[YY
][XX
] == 0.0
845 && inputrec
->deform
[YY
][YY
] == 0.0 && inputrec
->deform
[YY
][ZZ
] == 0.0
846 && inputrec
->deform
[ZZ
][XX
] == 0.0 && inputrec
->deform
[ZZ
][YY
] == 0.0
847 && inputrec
->deform
[ZZ
][ZZ
] == 0.0,
848 "Deformation is not supported by the modular simulator.");
851 && conditionalAssert(gmx_mtop_interaction_count(globalTopology
, IF_VSITE
) == 0,
852 "Virtual sites are not supported by the modular simulator.");
853 isInputCompatible
= isInputCompatible
854 && conditionalAssert(!inputrec
->bDoAwh
,
855 "AWH is not supported by the modular simulator.");
858 && conditionalAssert(gmx_mtop_ftype_count(globalTopology
, F_DISRES
) == 0,
859 "Distance restraints are not supported by the modular simulator.");
862 && conditionalAssert(
863 gmx_mtop_ftype_count(globalTopology
, F_ORIRES
) == 0,
864 "Orientation restraints are not supported by the modular simulator.");
867 && conditionalAssert(ms
== nullptr,
868 "Multi-sim are not supported by the modular simulator.");
871 && conditionalAssert(replExParams
.exchangeInterval
== 0,
872 "Replica exchange is not supported by the modular simulator.");
874 int numEnsembleRestraintSystems
;
877 numEnsembleRestraintSystems
= fcd
->disres
.nsystems
;
881 auto distantRestraintEnsembleEnvVar
= getenv("GMX_DISRE_ENSEMBLE_SIZE");
882 numEnsembleRestraintSystems
=
883 (ms
!= nullptr && distantRestraintEnsembleEnvVar
!= nullptr)
884 ? static_cast<int>(strtol(distantRestraintEnsembleEnvVar
, nullptr, 10))
889 && conditionalAssert(numEnsembleRestraintSystems
<= 1,
890 "Ensemble restraints are not supported by the modular simulator.");
893 && conditionalAssert(!doSimulatedAnnealing(inputrec
),
894 "Simulated annealing is not supported by the modular simulator.");
897 && conditionalAssert(!inputrec
->bSimTemp
,
898 "Simulated tempering is not supported by the modular simulator.");
899 isInputCompatible
= isInputCompatible
900 && conditionalAssert(!inputrec
->bExpanded
,
901 "Expanded ensemble simulations are not supported by "
902 "the modular simulator.");
905 && conditionalAssert(!doEssentialDynamics
,
906 "Essential dynamics is not supported by the modular simulator.");
907 isInputCompatible
= isInputCompatible
908 && conditionalAssert(inputrec
->eSwapCoords
== eswapNO
,
909 "Ion / water position swapping is not supported by "
910 "the modular simulator.");
913 && conditionalAssert(!inputrec
->bIMD
,
914 "Interactive MD is not supported by the modular simulator.");
917 && conditionalAssert(!doMembed
,
918 "Membrane embedding is not supported by the modular simulator.");
919 // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
922 && conditionalAssert(
923 getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
924 "Integration on the GPU is not supported by the modular simulator.");
925 // Modular simulator is centered around NS updates
926 // TODO: think how to handle nstlist == 0
927 isInputCompatible
= isInputCompatible
928 && conditionalAssert(inputrec
->nstlist
!= 0,
929 "Simulations without neighbor list update are not "
930 "supported by the modular simulator.");
931 isInputCompatible
= isInputCompatible
932 && conditionalAssert(!GMX_FAHCORE
,
933 "GMX_FAHCORE not supported by the modular simulator.");
935 return isInputCompatible
;
938 void ModularSimulator::checkInputForDisabledFunctionality()
940 isInputCompatible(true, inputrec
, doRerun
, *top_global
, ms
, replExParams
, fcd
,
941 opt2bSet("-ei", nfile
, fnm
), membed
!= nullptr);
942 if (observablesHistory
->edsamHistory
)
945 "The checkpoint is from a run with essential dynamics sampling, "
946 "but the current run did not specify the -ei option. "
947 "Either specify the -ei option to mdrun, or do not use this checkpoint file.");
951 SignallerCallbackPtr
ModularSimulator::SignalHelper::registerLastStepCallback()
953 return std::make_unique
<SignallerCallback
>(
954 [this](Step step
, Time gmx_unused time
) { this->lastStep_
= step
; });
957 SignallerCallbackPtr
ModularSimulator::SignalHelper::registerNSCallback()
959 return std::make_unique
<SignallerCallback
>(
960 [this](Step step
, Time gmx_unused time
) { this->nextNSStep_
= step
; });