Use more forward declarations to removed header dependencies
[alexxy/gromacs.git] / 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/math/vec.h"
54 #include "gromacs/mdlib/checkpointhandler.h"
55 #include "gromacs/mdlib/constr.h"
56 #include "gromacs/mdlib/energyoutput.h"
57 #include "gromacs/mdlib/mdatoms.h"
58 #include "gromacs/mdlib/resethandler.h"
59 #include "gromacs/mdlib/stat.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/mdtypes/state.h"
73 #include "gromacs/nbnxm/nbnxm.h"
74 #include "gromacs/timing/walltime_accounting.h"
75 #include "gromacs/topology/topology.h"
76 #include "gromacs/utility/cstringutil.h"
77 #include "gromacs/utility/fatalerror.h"
78
79 #include "compositesimulatorelement.h"
80 #include "computeglobalselement.h"
81 #include "constraintelement.h"
82 #include "energyelement.h"
83 #include "forceelement.h"
84 #include "freeenergyperturbationelement.h"
85 #include "parrinellorahmanbarostat.h"
86 #include "propagator.h"
87 #include "shellfcelement.h"
88 #include "signallers.h"
89 #include "statepropagatordata.h"
90 #include "trajectoryelement.h"
91 #include "vrescalethermostat.h"
92
93 namespace gmx
94 {
95 void ModularSimulator::run()
96 {
97     GMX_LOG(mdlog.info).asParagraph().appendText("Using the modular simulator.");
98     constructElementsAndSignallers();
99     simulatorSetup();
100     for (auto& signaller : signallerCallList_)
101     {
102         signaller->signallerSetup();
103     }
104     if (domDecHelper_)
105     {
106         domDecHelper_->setup();
107     }
108
109     for (auto& element : elementsOwnershipList_)
110     {
111         element->elementSetup();
112     }
113     if (pmeLoadBalanceHelper_)
114     {
115         // State must have been initialized so pmeLoadBalanceHelper_ gets a valid box
116         pmeLoadBalanceHelper_->setup();
117     }
118
119     while (step_ <= signalHelper_->lastStep_)
120     {
121         populateTaskQueue();
122
123         while (!taskQueue_.empty())
124         {
125             auto task = std::move(taskQueue_.front());
126             taskQueue_.pop();
127             // run function
128             (*task)();
129         }
130     }
131
132     for (auto& element : elementsOwnershipList_)
133     {
134         element->elementTeardown();
135     }
136     if (pmeLoadBalanceHelper_)
137     {
138         pmeLoadBalanceHelper_->teardown();
139     }
140     simulatorTeardown();
141 }
142
143 void ModularSimulator::simulatorSetup()
144 {
145     if (!mdrunOptions.writeConfout)
146     {
147         // This is on by default, and the main known use case for
148         // turning it off is for convenience in benchmarking, which is
149         // something that should not show up in the general user
150         // interface.
151         GMX_LOG(mdlog.info)
152                 .asParagraph()
153                 .appendText(
154                         "The -noconfout functionality is deprecated, and "
155                         "may be removed in a future version.");
156     }
157
158     if (MASTER(cr))
159     {
160         char        sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
161         std::string timeString;
162         fprintf(stderr, "starting mdrun '%s'\n", *(top_global->name));
163         if (inputrec->nsteps >= 0)
164         {
165             timeString = formatString("%8.1f", static_cast<double>(inputrec->init_step + inputrec->nsteps)
166                                                        * inputrec->delta_t);
167         }
168         else
169         {
170             timeString = "infinite";
171         }
172         if (inputrec->init_step > 0)
173         {
174             fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
175                     gmx_step_str(inputrec->init_step + inputrec->nsteps, sbuf), timeString.c_str(),
176                     gmx_step_str(inputrec->init_step, sbuf2), inputrec->init_step * inputrec->delta_t);
177         }
178         else
179         {
180             fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(inputrec->nsteps, sbuf),
181                     timeString.c_str());
182         }
183         fprintf(fplog, "\n");
184     }
185
186     walltime_accounting_start_time(walltime_accounting);
187     wallcycle_start(wcycle, ewcRUN);
188     print_start(fplog, cr, walltime_accounting, "mdrun");
189
190     step_ = inputrec->init_step;
191 }
192
193 void ModularSimulator::preStep(Step step, Time gmx_unused time, bool isNeighborSearchingStep)
194 {
195     if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) && step != signalHelper_->lastStep_)
196     {
197         /*
198          * Stop handler wants to stop after the current step, which was
199          * not known when building the current task queue. This happens
200          * e.g. when a stop is signalled by OS. We therefore want to purge
201          * the task queue now, and re-schedule this step as last step.
202          */
203         // clear task queue
204         std::queue<SimulatorRunFunctionPtr>().swap(taskQueue_);
205         // rewind step
206         step_ = step;
207         return;
208     }
209
210     resetHandler_->setSignal(walltime_accounting);
211     // This is a hack to avoid having to rewrite StopHandler to be a NeighborSearchSignaller
212     // and accept the step as input. Eventually, we want to do that, but currently this would
213     // require introducing NeighborSearchSignaller in the legacy do_md or a lot of code
214     // duplication.
215     stophandlerIsNSStep_    = isNeighborSearchingStep;
216     stophandlerCurrentStep_ = step;
217     stopHandler_->setSignal();
218
219     wallcycle_start(wcycle, ewcSTEP);
220 }
221
222 void ModularSimulator::postStep(Step step, Time gmx_unused time)
223 {
224     // Output stuff
225     if (MASTER(cr))
226     {
227         if (do_per_step(step, inputrec->nstlog))
228         {
229             if (fflush(fplog) != 0)
230             {
231                 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
232             }
233         }
234     }
235     const bool do_verbose = mdrunOptions.verbose
236                             && (step % mdrunOptions.verboseStepPrintInterval == 0
237                                 || step == inputrec->init_step || step == signalHelper_->lastStep_);
238     // Print the remaining wall clock time for the run
239     if (MASTER(cr) && (do_verbose || gmx_got_usr_signal())
240         && !(pmeLoadBalanceHelper_ && pmeLoadBalanceHelper_->pmePrinting()))
241     {
242         print_time(stderr, walltime_accounting, step, inputrec, cr);
243     }
244
245     double cycles = wallcycle_stop(wcycle, ewcSTEP);
246     if (DOMAINDECOMP(cr) && wcycle)
247     {
248         dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
249     }
250
251     resetHandler_->resetCounters(
252             step, step - inputrec->init_step, mdlog, fplog, cr, fr->nbv.get(), nrnb, fr->pmedata,
253             pmeLoadBalanceHelper_ ? pmeLoadBalanceHelper_->loadBalancingObject() : nullptr, wcycle,
254             walltime_accounting);
255 }
256
257 void ModularSimulator::simulatorTeardown()
258 {
259
260     // Stop measuring walltime
261     walltime_accounting_end_time(walltime_accounting);
262
263     if (!thisRankHasDuty(cr, DUTY_PME))
264     {
265         /* Tell the PME only node to finish */
266         gmx_pme_send_finish(cr);
267     }
268
269     walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
270 }
271
272 void ModularSimulator::populateTaskQueue()
273 {
274     auto registerRunFunction = std::make_unique<RegisterRunFunction>(
275             [this](SimulatorRunFunctionPtr ptr) { taskQueue_.push(std::move(ptr)); });
276
277     Time startTime = inputrec->init_t;
278     Time timeStep  = inputrec->delta_t;
279     Time time      = startTime + step_ * timeStep;
280
281     // Run an initial call to the signallers
282     for (auto& signaller : signallerCallList_)
283     {
284         signaller->signal(step_, time);
285     }
286
287     if (checkpointHelper_)
288     {
289         checkpointHelper_->run(step_, time);
290     }
291
292     if (pmeLoadBalanceHelper_)
293     {
294         pmeLoadBalanceHelper_->run(step_, time);
295     }
296     if (domDecHelper_)
297     {
298         domDecHelper_->run(step_, time);
299     }
300
301     do
302     {
303         // local variables for lambda capturing
304         const int  step     = step_;
305         const bool isNSStep = step == signalHelper_->nextNSStep_;
306
307         // register pre-step
308         (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
309                 [this, step, time, isNSStep]() { preStep(step, time, isNSStep); }));
310         // register elements for step
311         for (auto& element : elementCallList_)
312         {
313             element->scheduleTask(step_, time, registerRunFunction);
314         }
315         // register post-step
316         (*registerRunFunction)(
317                 std::make_unique<SimulatorRunFunction>([this, step, time]() { postStep(step, time); }));
318
319         // prepare next step
320         step_++;
321         time = startTime + step_ * timeStep;
322         for (auto& signaller : signallerCallList_)
323         {
324             signaller->signal(step_, time);
325         }
326     } while (step_ != signalHelper_->nextNSStep_ && step_ <= signalHelper_->lastStep_);
327 }
328
329 void ModularSimulator::constructElementsAndSignallers()
330 {
331     /* When restarting from a checkpoint, it can be appropriate to
332      * initialize ekind from quantities in the checkpoint. Otherwise,
333      * compute_globals must initialize ekind before the simulation
334      * starts/restarts. However, only the master rank knows what was
335      * found in the checkpoint file, so we have to communicate in
336      * order to coordinate the restart.
337      *
338      * TODO (modular) This should become obsolete when checkpoint reading
339      *      happens within the modular simulator framework: The energy
340      *      element should read its data from the checkpoint file pointer,
341      *      and signal to the compute globals element if it needs anything
342      *      reduced.
343      *
344      * TODO (legacy) Consider removing this communication if/when checkpoint
345      *      reading directly follows .tpr reading, because all ranks can
346      *      agree on hasReadEkinState at that time.
347      */
348     bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
349     if (PAR(cr))
350     {
351         gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
352     }
353     if (hasReadEkinState)
354     {
355         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
356     }
357
358     /*
359      * Build data structures
360      */
361     topologyHolder_ =
362             std::make_unique<TopologyHolder>(*top_global, cr, inputrec, fr, mdAtoms, constr, vsite);
363
364     std::unique_ptr<FreeEnergyPerturbationElement> freeEnergyPerturbationElement    = nullptr;
365     FreeEnergyPerturbationElement*                 freeEnergyPerturbationElementPtr = nullptr;
366     if (inputrec->efep != efepNO)
367     {
368         freeEnergyPerturbationElement =
369                 std::make_unique<FreeEnergyPerturbationElement>(fplog, inputrec, mdAtoms);
370         freeEnergyPerturbationElementPtr = freeEnergyPerturbationElement.get();
371     }
372
373     auto statePropagatorData = std::make_unique<StatePropagatorData>(
374             top_global->natoms, fplog, cr, state_global, inputrec->nstxout, inputrec->nstvout,
375             inputrec->nstfout, inputrec->nstxout_compressed, fr->nbv->useGpu(),
376             freeEnergyPerturbationElementPtr, topologyHolder_.get(), fr->bMolPBC,
377             mdrunOptions.writeConfout, opt2fn("-c", nfile, fnm), inputrec, mdAtoms->mdatoms());
378     auto statePropagatorDataPtr = compat::make_not_null(statePropagatorData.get());
379
380     auto energyElement = std::make_unique<EnergyElement>(
381             statePropagatorDataPtr, freeEnergyPerturbationElementPtr, top_global, inputrec, mdAtoms,
382             enerd, ekind, constr, fplog, fcd, mdModulesNotifier, MASTER(cr), observablesHistory,
383             startingBehavior);
384     auto energyElementPtr = compat::make_not_null(energyElement.get());
385
386     /*
387      * Build stop handler
388      */
389     const bool simulationsShareState = false;
390     stopHandler_                     = stopHandlerBuilder->getStopHandlerMD(
391             compat::not_null<SimulationSignal*>(&signals_[eglsSTOPCOND]), simulationsShareState,
392             MASTER(cr), inputrec->nstlist, mdrunOptions.reproducible, nstglobalcomm_,
393             mdrunOptions.maximumHoursToRun, inputrec->nstlist == 0, fplog, stophandlerCurrentStep_,
394             stophandlerIsNSStep_, walltime_accounting);
395
396     /*
397      * Create simulator builders
398      */
399     SignallerBuilder<NeighborSearchSignaller> neighborSearchSignallerBuilder;
400     SignallerBuilder<LastStepSignaller>       lastStepSignallerBuilder;
401     SignallerBuilder<LoggingSignaller>        loggingSignallerBuilder;
402     SignallerBuilder<EnergySignaller>         energySignallerBuilder;
403     TrajectoryElementBuilder                  trajectoryElementBuilder;
404
405     /*
406      * Register data structures to signallers
407      */
408     trajectoryElementBuilder.registerWriterClient(statePropagatorDataPtr);
409     trajectoryElementBuilder.registerSignallerClient(statePropagatorDataPtr);
410     lastStepSignallerBuilder.registerSignallerClient(statePropagatorDataPtr);
411
412     trajectoryElementBuilder.registerWriterClient(energyElementPtr);
413     trajectoryElementBuilder.registerSignallerClient(energyElementPtr);
414     energySignallerBuilder.registerSignallerClient(energyElementPtr);
415
416     // Register the simulator itself to the neighbor search / last step signaller
417     neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
418     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
419
420     /*
421      * Build integrator - this takes care of force calculation, propagation,
422      * constraining, and of the place the statePropagatorData and the energy element
423      * have a full timestep state.
424      */
425     // TODO: Make a CheckpointHelperBuilder
426     std::vector<ICheckpointHelperClient*> checkpointClients = { statePropagatorDataPtr, energyElementPtr,
427                                                                 freeEnergyPerturbationElementPtr };
428     CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback = nullptr;
429     auto                               integrator =
430             buildIntegrator(&neighborSearchSignallerBuilder, &energySignallerBuilder,
431                             &loggingSignallerBuilder, &trajectoryElementBuilder, &checkpointClients,
432                             &checkBondedInteractionsCallback, statePropagatorDataPtr,
433                             energyElementPtr, freeEnergyPerturbationElementPtr, hasReadEkinState);
434
435     /*
436      * Build infrastructure elements
437      */
438
439     if (PmeLoadBalanceHelper::doPmeLoadBalancing(mdrunOptions, inputrec, fr))
440     {
441         pmeLoadBalanceHelper_ = std::make_unique<PmeLoadBalanceHelper>(
442                 mdrunOptions.verbose, statePropagatorDataPtr, fplog, cr, mdlog, inputrec, wcycle, fr);
443         neighborSearchSignallerBuilder.registerSignallerClient(
444                 compat::make_not_null(pmeLoadBalanceHelper_.get()));
445     }
446
447     if (DOMAINDECOMP(cr))
448     {
449         GMX_ASSERT(checkBondedInteractionsCallback,
450                    "Domain decomposition needs a callback for check the number of bonded "
451                    "interactions.");
452         domDecHelper_ = std::make_unique<DomDecHelper>(
453                 mdrunOptions.verbose, mdrunOptions.verboseStepPrintInterval, statePropagatorDataPtr,
454                 topologyHolder_.get(), std::move(checkBondedInteractionsCallback), nstglobalcomm_, fplog,
455                 cr, mdlog, constr, inputrec, mdAtoms, nrnb, wcycle, fr, vsite, imdSession, pull_work);
456         neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(domDecHelper_.get()));
457     }
458
459     const bool simulationsShareResetCounters = false;
460     resetHandler_                            = std::make_unique<ResetHandler>(
461             compat::make_not_null<SimulationSignal*>(&signals_[eglsRESETCOUNTERS]),
462             simulationsShareResetCounters, inputrec->nsteps, MASTER(cr),
463             mdrunOptions.timingOptions.resetHalfway, mdrunOptions.maximumHoursToRun, mdlog, wcycle,
464             walltime_accounting);
465
466     /*
467      * Build signaller list
468      *
469      * Note that as signallers depend on each others, the order of calling the signallers
470      * matters. It is the responsibility of this builder to ensure that the order is
471      * maintained.
472      */
473     auto energySignaller = energySignallerBuilder.build(
474             inputrec->nstcalcenergy, inputrec->fepvals->nstdhdl, inputrec->nstpcouple);
475     trajectoryElementBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
476     loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
477     auto trajectoryElement = trajectoryElementBuilder.build(
478             fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec,
479             top_global, oenv, wcycle, startingBehavior);
480     loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
481
482     // Add checkpoint helper here since we need a pointer to the trajectory element and
483     // need to register it with the lastStepSignallerBuilder
484     auto checkpointHandler = std::make_unique<CheckpointHandler>(
485             compat::make_not_null<SimulationSignal*>(&signals_[eglsCHKPT]), simulationsShareState,
486             inputrec->nstlist == 0, MASTER(cr), mdrunOptions.writeConfout,
487             mdrunOptions.checkpointOptions.period);
488     checkpointHelper_ = std::make_unique<CheckpointHelper>(
489             std::move(checkpointClients), std::move(checkpointHandler), inputrec->init_step,
490             trajectoryElement.get(), top_global->natoms, fplog, cr, observablesHistory,
491             walltime_accounting, state_global, mdrunOptions.writeConfout);
492     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(checkpointHelper_.get()));
493
494     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
495     auto loggingSignaller =
496             loggingSignallerBuilder.build(inputrec->nstlog, inputrec->init_step, inputrec->init_t);
497     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(loggingSignaller.get()));
498     auto lastStepSignaller =
499             lastStepSignallerBuilder.build(inputrec->nsteps, inputrec->init_step, stopHandler_.get());
500     neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(lastStepSignaller.get()));
501     auto neighborSearchSignaller = neighborSearchSignallerBuilder.build(
502             inputrec->nstlist, inputrec->init_step, inputrec->init_t);
503
504     addToCallListAndMove(std::move(neighborSearchSignaller), signallerCallList_, signallersOwnershipList_);
505     addToCallListAndMove(std::move(lastStepSignaller), signallerCallList_, signallersOwnershipList_);
506     addToCallListAndMove(std::move(loggingSignaller), signallerCallList_, signallersOwnershipList_);
507     addToCallList(trajectoryElement, signallerCallList_);
508     addToCallListAndMove(std::move(energySignaller), signallerCallList_, signallersOwnershipList_);
509
510     /*
511      * Build the element list
512      *
513      * This is the actual sequence of (non-infrastructure) elements to be run.
514      * For NVE, the trajectory element is used outside of the integrator
515      * (composite) element, as well as the checkpoint helper. The checkpoint
516      * helper should be on top of the loop, and is only part of the simulator
517      * call list to be able to react to the last step being signalled.
518      */
519     addToCallList(checkpointHelper_, elementCallList_);
520     if (freeEnergyPerturbationElement)
521     {
522         addToCallListAndMove(std::move(freeEnergyPerturbationElement), elementCallList_,
523                              elementsOwnershipList_);
524     }
525     addToCallListAndMove(std::move(integrator), elementCallList_, elementsOwnershipList_);
526     addToCallListAndMove(std::move(trajectoryElement), elementCallList_, elementsOwnershipList_);
527     // for vv, we need to setup statePropagatorData after the compute
528     // globals so that we reset the right velocities
529     // TODO: Avoid this by getting rid of the need of resetting velocities in vv
530     elementsOwnershipList_.emplace_back(std::move(statePropagatorData));
531     elementsOwnershipList_.emplace_back(std::move(energyElement));
532 }
533
534 std::unique_ptr<ISimulatorElement>
535 ModularSimulator::buildForces(SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
536                               SignallerBuilder<EnergySignaller>*         energySignallerBuilder,
537                               StatePropagatorData*                       statePropagatorDataPtr,
538                               EnergyElement*                             energyElementPtr,
539                               FreeEnergyPerturbationElement* freeEnergyPerturbationElement)
540 {
541     const bool isVerbose    = mdrunOptions.verbose;
542     const bool isDynamicBox = inputrecDynamicBox(inputrec);
543     // Check for polarizable models and flexible constraints
544     if (ShellFCElement::doShellsOrFlexConstraints(topologyHolder_->globalTopology(),
545                                                   constr ? constr->numFlexibleConstraints() : 0))
546     {
547         auto shellFCElement = std::make_unique<ShellFCElement>(
548                 statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElement, isVerbose,
549                 isDynamicBox, fplog, cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, runScheduleWork, vsite,
550                 imdSession, pull_work, constr, &topologyHolder_->globalTopology(), enforcedRotation);
551         topologyHolder_->registerClient(shellFCElement.get());
552         neighborSearchSignallerBuilder->registerSignallerClient(
553                 compat::make_not_null(shellFCElement.get()));
554         energySignallerBuilder->registerSignallerClient(compat::make_not_null(shellFCElement.get()));
555
556         // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
557         return std::move(shellFCElement);
558     }
559     else
560     {
561         auto forceElement = std::make_unique<ForceElement>(
562                 statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElement,
563                 isDynamicBox, fplog, cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, runScheduleWork,
564                 vsite, imdSession, pull_work, enforcedRotation);
565         topologyHolder_->registerClient(forceElement.get());
566         neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
567         energySignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
568
569         // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
570         return std::move(forceElement);
571     }
572 }
573
574 std::unique_ptr<ISimulatorElement> ModularSimulator::buildIntegrator(
575         SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
576         SignallerBuilder<EnergySignaller>*         energySignallerBuilder,
577         SignallerBuilder<LoggingSignaller>*        loggingSignallerBuilder,
578         TrajectoryElementBuilder*                  trajectoryElementBuilder,
579         std::vector<ICheckpointHelperClient*>*     checkpointClients,
580         CheckBondedInteractionsCallbackPtr*        checkBondedInteractionsCallback,
581         compat::not_null<StatePropagatorData*>     statePropagatorDataPtr,
582         compat::not_null<EnergyElement*>           energyElementPtr,
583         FreeEnergyPerturbationElement*             freeEnergyPerturbationElementPtr,
584         bool                                       hasReadEkinState)
585 {
586     auto forceElement =
587             buildForces(neighborSearchSignallerBuilder, energySignallerBuilder,
588                         statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr);
589
590     // list of elements owned by the simulator composite object
591     std::vector<std::unique_ptr<ISimulatorElement>> elementsOwnershipList;
592     // call list of the simulator composite object
593     std::vector<compat::not_null<ISimulatorElement*>> elementCallList;
594
595     std::function<void()> needToCheckNumberOfBondedInteractions;
596     if (inputrec->eI == eiMD)
597     {
598         auto computeGlobalsElement =
599                 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>(
600                         statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
601                         &signals_, nstglobalcomm_, fplog, mdlog, cr, inputrec, mdAtoms, nrnb,
602                         wcycle, fr, &topologyHolder_->globalTopology(), constr, hasReadEkinState);
603         topologyHolder_->registerClient(computeGlobalsElement.get());
604         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
605         trajectoryElementBuilder->registerSignallerClient(
606                 compat::make_not_null(computeGlobalsElement.get()));
607
608         *checkBondedInteractionsCallback =
609                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
610
611         auto propagator = std::make_unique<Propagator<IntegrationStep::LeapFrog>>(
612                 inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
613
614         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
615         addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
616         if (inputrec->etc == etcVRESCALE)
617         {
618             // TODO: With increased complexity of the propagator, this will need further development,
619             //       e.g. using propagators templated for velocity propagation policies and a builder
620             propagator->setNumVelocityScalingVariables(inputrec->opts.ngtc);
621             auto thermostat = std::make_unique<VRescaleThermostat>(
622                     inputrec->nsttcouple, -1, false, inputrec->ld_seed, inputrec->opts.ngtc,
623                     inputrec->delta_t * inputrec->nsttcouple, inputrec->opts.ref_t, inputrec->opts.tau_t,
624                     inputrec->opts.nrdf, energyElementPtr, propagator->viewOnVelocityScaling(),
625                     propagator->velocityScalingCallback(), state_global, cr, inputrec->bContinuation);
626             checkpointClients->emplace_back(thermostat.get());
627             energyElementPtr->setVRescaleThermostat(thermostat.get());
628             addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
629         }
630
631         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
632         if (inputrec->epc == epcPARRINELLORAHMAN)
633         {
634             // Building the PR barostat here since it needs access to the propagator
635             // and we want to be able to move the propagator object
636             prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
637                     inputrec->nstpcouple, -1, inputrec->delta_t * inputrec->nstpcouple,
638                     inputrec->init_step, propagator->viewOnPRScalingMatrix(),
639                     propagator->prScalingCallback(), statePropagatorDataPtr, energyElementPtr,
640                     fplog, inputrec, mdAtoms, state_global, cr, inputrec->bContinuation);
641             energyElementPtr->setParrinelloRahamnBarostat(prBarostat.get());
642             checkpointClients->emplace_back(prBarostat.get());
643         }
644         addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
645         if (constr)
646         {
647             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
648                     constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
649                     MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
650             auto constraintElementPtr = compat::make_not_null(constraintElement.get());
651             energySignallerBuilder->registerSignallerClient(constraintElementPtr);
652             trajectoryElementBuilder->registerSignallerClient(constraintElementPtr);
653             loggingSignallerBuilder->registerSignallerClient(constraintElementPtr);
654
655             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
656         }
657
658         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
659         addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
660         if (prBarostat)
661         {
662             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
663         }
664     }
665     else if (inputrec->eI == eiVV)
666     {
667         auto computeGlobalsElement =
668                 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>(
669                         statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
670                         &signals_, nstglobalcomm_, fplog, mdlog, cr, inputrec, mdAtoms, nrnb,
671                         wcycle, fr, &topologyHolder_->globalTopology(), constr, hasReadEkinState);
672         topologyHolder_->registerClient(computeGlobalsElement.get());
673         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
674         trajectoryElementBuilder->registerSignallerClient(
675                 compat::make_not_null(computeGlobalsElement.get()));
676
677         *checkBondedInteractionsCallback =
678                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
679
680         auto propagatorVelocities = std::make_unique<Propagator<IntegrationStep::VelocitiesOnly>>(
681                 inputrec->delta_t * 0.5, statePropagatorDataPtr, mdAtoms, wcycle);
682         auto propagatorVelocitiesAndPositions =
683                 std::make_unique<Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>>(
684                         inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
685
686         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
687
688         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
689         if (inputrec->epc == epcPARRINELLORAHMAN)
690         {
691             // Building the PR barostat here since it needs access to the propagator
692             // and we want to be able to move the propagator object
693             prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
694                     inputrec->nstpcouple, -1, inputrec->delta_t * inputrec->nstpcouple,
695                     inputrec->init_step, propagatorVelocities->viewOnPRScalingMatrix(),
696                     propagatorVelocities->prScalingCallback(), statePropagatorDataPtr, energyElementPtr,
697                     fplog, inputrec, mdAtoms, state_global, cr, inputrec->bContinuation);
698             energyElementPtr->setParrinelloRahamnBarostat(prBarostat.get());
699             checkpointClients->emplace_back(prBarostat.get());
700         }
701         addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
702         if (constr)
703         {
704             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Velocities>>(
705                     constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
706                     MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
707             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
708             trajectoryElementBuilder->registerSignallerClient(
709                     compat::make_not_null(constraintElement.get()));
710             loggingSignallerBuilder->registerSignallerClient(
711                     compat::make_not_null(constraintElement.get()));
712
713             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
714         }
715         addToCallList(compat::make_not_null(computeGlobalsElement.get()), elementCallList);
716         addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
717         if (inputrec->etc == etcVRESCALE)
718         {
719             // TODO: With increased complexity of the propagator, this will need further development,
720             //       e.g. using propagators templated for velocity propagation policies and a builder
721             propagatorVelocitiesAndPositions->setNumVelocityScalingVariables(inputrec->opts.ngtc);
722             auto thermostat = std::make_unique<VRescaleThermostat>(
723                     inputrec->nsttcouple, 0, true, inputrec->ld_seed, inputrec->opts.ngtc,
724                     inputrec->delta_t * inputrec->nsttcouple, inputrec->opts.ref_t,
725                     inputrec->opts.tau_t, inputrec->opts.nrdf, energyElementPtr,
726                     propagatorVelocitiesAndPositions->viewOnVelocityScaling(),
727                     propagatorVelocitiesAndPositions->velocityScalingCallback(), state_global, cr,
728                     inputrec->bContinuation);
729             checkpointClients->emplace_back(thermostat.get());
730             energyElementPtr->setVRescaleThermostat(thermostat.get());
731             addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
732         }
733         addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList,
734                              elementsOwnershipList);
735         if (constr)
736         {
737             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
738                     constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
739                     MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
740             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
741             trajectoryElementBuilder->registerSignallerClient(
742                     compat::make_not_null(constraintElement.get()));
743             loggingSignallerBuilder->registerSignallerClient(
744                     compat::make_not_null(constraintElement.get()));
745
746             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
747         }
748         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
749         addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
750         if (prBarostat)
751         {
752             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
753         }
754     }
755     else
756     {
757         gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
758     }
759
760     auto integrator = std::make_unique<CompositeSimulatorElement>(std::move(elementCallList),
761                                                                   std::move(elementsOwnershipList));
762     // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
763     return std::move(integrator);
764 }
765
766 bool ModularSimulator::isInputCompatible(bool                             exitOnFailure,
767                                          const t_inputrec*                inputrec,
768                                          bool                             doRerun,
769                                          const gmx_vsite_t*               vsite,
770                                          const gmx_multisim_t*            ms,
771                                          const ReplicaExchangeParameters& replExParams,
772                                          const t_fcdata*                  fcd,
773                                          int                              nfile,
774                                          const t_filenm*                  fnm,
775                                          ObservablesHistory*              observablesHistory,
776                                          const gmx_membed_t*              membed)
777 {
778     auto conditionalAssert = [exitOnFailure](bool condition, const char* message) {
779         if (exitOnFailure)
780         {
781             GMX_RELEASE_ASSERT(condition, message);
782         }
783         return condition;
784     };
785
786     bool isInputCompatible = true;
787
788     // GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
789     // such as the leap-frog integrator
790     const auto modularSimulatorExplicitlyTurnedOn = (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
791     // GMX_USE_MODULAR_SIMULATOR allows to use disable modular simulator for all uses,
792     // including the velocity-verlet integrator used by default
793     const auto modularSimulatorExplicitlyTurnedOff = (getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
794
795     GMX_RELEASE_ASSERT(
796             !(modularSimulatorExplicitlyTurnedOn && modularSimulatorExplicitlyTurnedOff),
797             "Cannot have both GMX_USE_MODULAR_SIMULATOR=ON and GMX_DISABLE_MODULAR_SIMULATOR=ON. "
798             "Unset one of the two environment variables to explicitly chose which simulator to "
799             "use, "
800             "or unset both to recover default behavior.");
801
802     GMX_RELEASE_ASSERT(
803             !(modularSimulatorExplicitlyTurnedOff && inputrec->eI == eiVV
804               && inputrec->epc == epcPARRINELLORAHMAN),
805             "Cannot use a Parrinello-Rahman barostat with md-vv and "
806             "GMX_DISABLE_MODULAR_SIMULATOR=ON, "
807             "as the Parrinello-Rahman barostat is not implemented in the legacy simulator. Unset "
808             "GMX_DISABLE_MODULAR_SIMULATOR or use a different pressure control algorithm.");
809
810     isInputCompatible =
811             isInputCompatible
812             && conditionalAssert(
813                        inputrec->eI == eiMD || inputrec->eI == eiVV,
814                        "Only integrators md and md-vv are supported by the modular simulator.");
815     isInputCompatible = isInputCompatible
816                         && conditionalAssert(inputrec->eI != eiMD || modularSimulatorExplicitlyTurnedOn,
817                                              "Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular "
818                                              "simulator with integrator md.");
819     isInputCompatible =
820             isInputCompatible
821             && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
822     isInputCompatible =
823             isInputCompatible
824             && conditionalAssert(
825                        inputrec->etc == etcNO || inputrec->etc == etcVRESCALE,
826                        "Only v-rescale thermostat is supported by the modular simulator.");
827     isInputCompatible =
828             isInputCompatible
829             && conditionalAssert(
830                        inputrec->epc == epcNO || inputrec->epc == epcPARRINELLORAHMAN,
831                        "Only Parrinello-Rahman barostat is supported by the modular simulator.");
832     isInputCompatible =
833             isInputCompatible
834             && conditionalAssert(
835                        !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec)
836                          || inputrecNvtTrotter(inputrec)),
837                        "Legacy Trotter decomposition is not supported by the modular simulator.");
838     isInputCompatible = isInputCompatible
839                         && conditionalAssert(inputrec->efep == efepNO || inputrec->efep == efepYES
840                                                      || inputrec->efep == efepSLOWGROWTH,
841                                              "Expanded ensemble free energy calculation is not "
842                                              "supported by the modular simulator.");
843     isInputCompatible = isInputCompatible
844                         && conditionalAssert(!inputrec->bPull,
845                                              "Pulling is not supported by the modular simulator.");
846     isInputCompatible =
847             isInputCompatible
848             && conditionalAssert(inputrec->opts.ngacc == 1 && inputrec->opts.acc[0][XX] == 0.0
849                                          && inputrec->opts.acc[0][YY] == 0.0
850                                          && inputrec->opts.acc[0][ZZ] == 0.0 && inputrec->cos_accel == 0.0,
851                                  "Acceleration is not supported by the modular simulator.");
852     isInputCompatible =
853             isInputCompatible
854             && conditionalAssert(inputrec->opts.ngfrz == 1 && inputrec->opts.nFreeze[0][XX] == 0
855                                          && inputrec->opts.nFreeze[0][YY] == 0
856                                          && inputrec->opts.nFreeze[0][ZZ] == 0,
857                                  "Freeze groups are not supported by the modular simulator.");
858     isInputCompatible =
859             isInputCompatible
860             && conditionalAssert(
861                        inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
862                                && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
863                                && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
864                                && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
865                                && inputrec->deform[ZZ][ZZ] == 0.0,
866                        "Deformation is not supported by the modular simulator.");
867     isInputCompatible =
868             isInputCompatible
869             && conditionalAssert(vsite == nullptr,
870                                  "Virtual sites are not supported by the modular simulator.");
871     isInputCompatible = isInputCompatible
872                         && conditionalAssert(!inputrec->bDoAwh,
873                                              "AWH is not supported by the modular simulator.");
874     isInputCompatible =
875             isInputCompatible
876             && conditionalAssert(ms == nullptr,
877                                  "Multi-sim are not supported by the modular simulator.");
878     isInputCompatible =
879             isInputCompatible
880             && conditionalAssert(replExParams.exchangeInterval == 0,
881                                  "Replica exchange is not supported by the modular simulator.");
882     isInputCompatible =
883             isInputCompatible
884             && conditionalAssert(fcd->disres.nsystems <= 1,
885                                  "Ensemble restraints are not supported by the modular simulator.");
886     isInputCompatible =
887             isInputCompatible
888             && conditionalAssert(!doSimulatedAnnealing(inputrec),
889                                  "Simulated annealing is not supported by the modular simulator.");
890     isInputCompatible =
891             isInputCompatible
892             && conditionalAssert(!inputrec->bSimTemp,
893                                  "Simulated tempering is not supported by the modular simulator.");
894     isInputCompatible = isInputCompatible
895                         && conditionalAssert(!inputrec->bExpanded,
896                                              "Expanded ensemble simulations are not supported by "
897                                              "the modular simulator.");
898     isInputCompatible =
899             isInputCompatible
900             && conditionalAssert(
901                        !(opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr),
902                        "Essential dynamics is not supported by the modular simulator.");
903     isInputCompatible = isInputCompatible
904                         && conditionalAssert(inputrec->eSwapCoords == eswapNO,
905                                              "Ion / water position swapping is not supported by "
906                                              "the modular simulator.");
907     isInputCompatible =
908             isInputCompatible
909             && conditionalAssert(!inputrec->bIMD,
910                                  "Interactive MD is not supported by the modular simulator.");
911     isInputCompatible =
912             isInputCompatible
913             && conditionalAssert(membed == nullptr,
914                                  "Membrane embedding is not supported by the modular simulator.");
915     // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
916     isInputCompatible =
917             isInputCompatible
918             && conditionalAssert(
919                        getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
920                        "Integration on the GPU is not supported by the modular simulator.");
921     // Modular simulator is centered around NS updates
922     // TODO: think how to handle nstlist == 0
923     isInputCompatible = isInputCompatible
924                         && conditionalAssert(inputrec->nstlist != 0,
925                                              "Simulations without neighbor list update are not "
926                                              "supported by the modular simulator.");
927     isInputCompatible = isInputCompatible
928                         && conditionalAssert(!GMX_FAHCORE,
929                                              "GMX_FAHCORE not supported by the modular simulator.");
930
931     return isInputCompatible;
932 }
933
934 void ModularSimulator::checkInputForDisabledFunctionality()
935 {
936     isInputCompatible(true, inputrec, doRerun, vsite, ms, replExParams, fcd, nfile, fnm,
937                       observablesHistory, membed);
938 }
939
940 SignallerCallbackPtr ModularSimulator::SignalHelper::registerLastStepCallback()
941 {
942     return std::make_unique<SignallerCallback>(
943             [this](Step step, Time gmx_unused time) { this->lastStep_ = step; });
944 }
945
946 SignallerCallbackPtr ModularSimulator::SignalHelper::registerNSCallback()
947 {
948     return std::make_unique<SignallerCallback>(
949             [this](Step step, Time gmx_unused time) { this->nextNSStep_ = step; });
950 }
951 } // namespace gmx