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 global reduction element for the modular simulator
38 * \author Pascal Merz <pascal.merz@me.com>
39 * \ingroup module_modularsimulator
44 #include "computeglobalselement.h"
46 #include "gromacs/domdec/partition.h"
47 #include "gromacs/gmxlib/network.h"
48 #include "gromacs/gmxlib/nrnb.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdlib/md_support.h"
51 #include "gromacs/mdlib/mdatoms.h"
52 #include "gromacs/mdlib/stat.h"
53 #include "gromacs/mdlib/update.h"
54 #include "gromacs/mdtypes/commrec.h"
55 #include "gromacs/mdtypes/group.h"
56 #include "gromacs/mdtypes/inputrec.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/mdtypes/mdatom.h"
59 #include "gromacs/topology/topology.h"
61 #include "freeenergyperturbationdata.h"
62 #include "modularsimulator.h"
63 #include "simulatoralgorithm.h"
67 template<ComputeGlobalsAlgorithm algorithm>
68 ComputeGlobalsElement<algorithm>::ComputeGlobalsElement(StatePropagatorData* statePropagatorData,
69 EnergyData* energyData,
70 FreeEnergyPerturbationData* freeEnergyPerturbationData,
71 SimulationSignals* signals,
74 const MDLogger& mdlog,
76 const t_inputrec* inputrec,
77 const MDAtoms* mdAtoms,
79 gmx_wallcycle* wcycle,
81 const gmx_mtop_t* global_top,
82 Constraints* constr) :
83 energyReductionStep_(-1),
84 virialReductionStep_(-1),
85 vvSchedulingStep_(-1),
86 doStopCM_(inputrec->comm_mode != ecmNO),
87 nstcomm_(inputrec->nstcomm),
88 nstglobalcomm_(nstglobalcomm),
89 lastStep_(inputrec->nsteps + inputrec->init_step),
90 initStep_(inputrec->init_step),
91 nullSignaller_(std::make_unique<SimulationSignaller>(nullptr, nullptr, nullptr, false, false)),
92 totalNumberOfBondedInteractions_(0),
93 shouldCheckNumberOfBondedInteractions_(false),
94 statePropagatorData_(statePropagatorData),
95 energyData_(energyData),
96 localTopology_(nullptr),
97 freeEnergyPerturbationData_(freeEnergyPerturbationData),
98 vcm_(global_top->groups, *inputrec),
104 top_global_(global_top),
111 reportComRemovalInfo(fplog, vcm_);
112 gstat_ = global_stat_init(inputrec_);
115 template<ComputeGlobalsAlgorithm algorithm>
116 ComputeGlobalsElement<algorithm>::~ComputeGlobalsElement()
118 global_stat_destroy(gstat_);
121 template<ComputeGlobalsAlgorithm algorithm>
122 void ComputeGlobalsElement<algorithm>::elementSetup()
124 GMX_ASSERT(localTopology_, "Setup called before local topology was set.");
126 if (doStopCM_ && !inputrec_->bContinuation)
128 // To minimize communication, compute_globals computes the COM velocity
129 // and the kinetic energy for the velocities without COM motion removed.
130 // Thus to get the kinetic energy without the COM contribution, we need
131 // to call compute_globals twice.
133 compute(-1, CGLO_GSTAT | CGLO_STOPCM, nullSignaller_.get(), false, true);
135 auto v = statePropagatorData_->velocitiesView();
136 // At initialization, do not pass x with acceleration-correction mode
137 // to avoid (incorrect) correction of the initial coordinates.
138 auto x = vcm_.mode == ecmLINEAR_ACCELERATION_CORRECTION ? ArrayRefWithPadding<RVec>()
139 : statePropagatorData_->positionsView();
140 process_and_stopcm_grp(fplog_, &vcm_, *mdAtoms_->mdatoms(), x.unpaddedArrayRef(),
141 v.unpaddedArrayRef());
142 inc_nrnb(nrnb_, eNR_STOPCM, mdAtoms_->mdatoms()->homenr);
145 unsigned int cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
146 | (energyData_->hasReadEkinFromCheckpoint() ? CGLO_READEKIN : 0));
148 if (algorithm == ComputeGlobalsAlgorithm::VelocityVerlet)
150 cglo_flags |= CGLO_PRESSURE | CGLO_CONSTRAINT;
153 compute(-1, cglo_flags, nullSignaller_.get(), false, true);
155 // Calculate the initial half step temperature, and save the ekinh_old
156 for (int i = 0; (i < inputrec_->opts.ngtc); i++)
158 copy_mat(energyData_->ekindata()->tcstat[i].ekinh, energyData_->ekindata()->tcstat[i].ekinh_old);
162 template<ComputeGlobalsAlgorithm algorithm>
163 void ComputeGlobalsElement<algorithm>::scheduleTask(Step step,
164 Time gmx_unused time,
165 const RegisterRunFunction& registerRunFunction)
167 const bool needComReduction = doStopCM_ && do_per_step(step, nstcomm_);
168 const bool needGlobalReduction = step == energyReductionStep_ || step == virialReductionStep_
169 || needComReduction || do_per_step(step, nstglobalcomm_);
171 // TODO: CGLO_GSTAT is only used for needToSumEkinhOld_, i.e. to signal that we do or do not
172 // sum the previous kinetic energy. We should simplify / clarify this.
174 if (algorithm == ComputeGlobalsAlgorithm::LeapFrog)
176 // With Leap-Frog we can skip compute_globals at
177 // non-communication steps, but we need to calculate
178 // the kinetic energy one step before communication.
180 // With leap-frog we also need to compute the half-step
181 // kinetic energy at the step before we need to write
182 // the full-step kinetic energy
183 const bool needEkinAtNextStep = (do_per_step(step + 1, nstglobalcomm_) || step + 1 == lastStep_);
185 if (!needGlobalReduction && !needEkinAtNextStep)
190 const bool doEnergy = step == energyReductionStep_;
191 int flags = (needGlobalReduction ? CGLO_GSTAT : 0) | (doEnergy ? CGLO_ENERGY : 0)
192 | (needComReduction ? CGLO_STOPCM : 0) | CGLO_TEMPERATURE | CGLO_PRESSURE
195 // Since we're already communicating at this step, we
196 // can propagate intra-simulation signals. Note that
197 // check_nstglobalcomm has the responsibility for
198 // choosing the value of nstglobalcomm which satisfies
199 // the need of the different signallers.
200 const bool doIntraSimSignal = true;
201 // Disable functionality
202 const bool doInterSimSignal = false;
204 // Make signaller to signal stop / reset / checkpointing signals
205 auto signaller = std::make_shared<SimulationSignaller>(signals_, cr_, nullptr,
206 doInterSimSignal, doIntraSimSignal);
208 registerRunFunction([this, step, flags, signaller = std::move(signaller)]() {
209 compute(step, flags, signaller.get(), true);
212 else if (algorithm == ComputeGlobalsAlgorithm::VelocityVerlet)
214 // For VV, we schedule two calls to compute globals per step.
215 if (step != vvSchedulingStep_)
217 // This is the first scheduling call for this step (positions & velocities at full time
218 // step) Set this as the current scheduling step
219 vvSchedulingStep_ = step;
221 // For vv, the state at the beginning of the step is positions at time t, velocities at time t - dt/2
222 // The first velocity propagation (+dt/2) therefore actually corresponds to the previous step.
223 // So we need information from the last step in the first half of the integration
224 if (!needGlobalReduction && !do_per_step(step - 1, nstglobalcomm_))
229 const bool doTemperature = step != initStep_ || inputrec_->bContinuation;
230 const bool doEnergy = step == energyReductionStep_;
232 int flags = (needGlobalReduction ? CGLO_GSTAT : 0) | (doEnergy ? CGLO_ENERGY : 0)
233 | (doTemperature ? CGLO_TEMPERATURE : 0) | CGLO_PRESSURE | CGLO_CONSTRAINT
234 | (needComReduction ? CGLO_STOPCM : 0) | CGLO_SCALEEKIN;
237 [this, step, flags]() { compute(step, flags, nullSignaller_.get(), false); });
241 // second call to compute_globals for this step
242 // Reset the scheduling step to avoid confusion if scheduling needs
243 // to be repeated (in case of unexpected simulation termination)
244 vvSchedulingStep_ = -1;
246 if (!needGlobalReduction)
250 int flags = CGLO_GSTAT | CGLO_CONSTRAINT;
252 // Since we're already communicating at this step, we
253 // can propagate intra-simulation signals. Note that
254 // check_nstglobalcomm has the responsibility for
255 // choosing the value of nstglobalcomm which satisfies
256 // the need of the different signallers.
257 const bool doIntraSimSignal = true;
258 // Disable functionality
259 const bool doInterSimSignal = false;
261 auto signaller = std::make_shared<SimulationSignaller>(
262 signals_, cr_, nullptr, doInterSimSignal, doIntraSimSignal);
264 registerRunFunction([this, step, flags, signaller = std::move(signaller)]() {
265 compute(step, flags, signaller.get(), true);
271 template<ComputeGlobalsAlgorithm algorithm>
272 void ComputeGlobalsElement<algorithm>::compute(gmx::Step step,
274 SimulationSignaller* signaller,
278 auto x = statePropagatorData_->positionsView().unpaddedArrayRef();
279 auto v = statePropagatorData_->velocitiesView().unpaddedArrayRef();
280 auto box = statePropagatorData_->constBox();
281 auto lastbox = useLastBox ? statePropagatorData_->constPreviousBox()
282 : statePropagatorData_->constBox();
285 gstat_, cr_, inputrec_, fr_, energyData_->ekindata(), x, v, box, mdAtoms_->mdatoms(),
286 nrnb_, &vcm_, step != -1 ? wcycle_ : nullptr, energyData_->enerdata(),
287 energyData_->forceVirial(step), energyData_->constraintVirial(step),
288 energyData_->totalVirial(step), energyData_->pressure(step), constr_, signaller,
289 lastbox, &totalNumberOfBondedInteractions_, energyData_->needToSumEkinhOld(),
290 flags | (shouldCheckNumberOfBondedInteractions_ ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
291 checkNumberOfBondedInteractions(mdlog_, cr_, totalNumberOfBondedInteractions_, top_global_,
292 localTopology_, x, box, &shouldCheckNumberOfBondedInteractions_);
293 if (flags & CGLO_STOPCM && !isInit)
295 process_and_stopcm_grp(fplog_, &vcm_, *mdAtoms_->mdatoms(), x, v);
296 inc_nrnb(nrnb_, eNR_STOPCM, mdAtoms_->mdatoms()->homenr);
300 template<ComputeGlobalsAlgorithm algorithm>
301 CheckBondedInteractionsCallback ComputeGlobalsElement<algorithm>::getCheckNumberOfBondedInteractionsCallback()
303 return [this]() { needToCheckNumberOfBondedInteractions(); };
306 template<ComputeGlobalsAlgorithm algorithm>
307 void ComputeGlobalsElement<algorithm>::needToCheckNumberOfBondedInteractions()
309 shouldCheckNumberOfBondedInteractions_ = true;
312 template<ComputeGlobalsAlgorithm algorithm>
313 void ComputeGlobalsElement<algorithm>::setTopology(const gmx_localtop_t* top)
315 localTopology_ = top;
318 template<ComputeGlobalsAlgorithm algorithm>
319 std::optional<SignallerCallback> ComputeGlobalsElement<algorithm>::registerEnergyCallback(EnergySignallerEvent event)
321 if (event == EnergySignallerEvent::EnergyCalculationStep)
323 return [this](Step step, Time /*unused*/) { energyReductionStep_ = step; };
325 if (event == EnergySignallerEvent::VirialCalculationStep)
327 return [this](Step step, Time /*unused*/) { virialReductionStep_ = step; };
332 template<ComputeGlobalsAlgorithm algorithm>
333 std::optional<SignallerCallback>
334 ComputeGlobalsElement<algorithm>::registerTrajectorySignallerCallback(TrajectoryEvent event)
336 if (event == TrajectoryEvent::EnergyWritingStep)
338 return [this](Step step, Time /*unused*/) { energyReductionStep_ = step; };
343 //! Explicit template instantiation
345 template class ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>;
346 template class ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>;
350 ISimulatorElement* ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>::getElementPointerImpl(
351 LegacySimulatorData* legacySimulatorData,
352 ModularSimulatorAlgorithmBuilderHelper* builderHelper,
353 StatePropagatorData* statePropagatorData,
354 EnergyData* energyData,
355 FreeEnergyPerturbationData* freeEnergyPerturbationData,
356 GlobalCommunicationHelper* globalCommunicationHelper)
358 auto* element = builderHelper->storeElement(
359 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>(
360 statePropagatorData, energyData, freeEnergyPerturbationData,
361 globalCommunicationHelper->simulationSignals(),
362 globalCommunicationHelper->nstglobalcomm(), legacySimulatorData->fplog,
363 legacySimulatorData->mdlog, legacySimulatorData->cr,
364 legacySimulatorData->inputrec, legacySimulatorData->mdAtoms,
365 legacySimulatorData->nrnb, legacySimulatorData->wcycle, legacySimulatorData->fr,
366 legacySimulatorData->top_global, legacySimulatorData->constr));
368 // TODO: Remove this when DD can reduce bonded interactions independently (#3421)
369 auto* castedElement = static_cast<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>*>(element);
370 globalCommunicationHelper->setCheckBondedInteractionsCallback(
371 castedElement->getCheckNumberOfBondedInteractionsCallback());
377 ISimulatorElement* ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>::getElementPointerImpl(
378 LegacySimulatorData* simulator,
379 ModularSimulatorAlgorithmBuilderHelper* builderHelper,
380 StatePropagatorData* statePropagatorData,
381 EnergyData* energyData,
382 FreeEnergyPerturbationData* freeEnergyPerturbationData,
383 GlobalCommunicationHelper* globalCommunicationHelper)
385 // We allow this element to be added multiple times to the call list, but we only want one
386 // actual element built
387 static const std::string key("vvComputeGlobalsElement");
389 const std::optional<std::any> cachedValue = builderHelper->getStoredValue(key);
393 return std::any_cast<ISimulatorElement*>(cachedValue.value());
397 ISimulatorElement* vvComputeGlobalsElement = builderHelper->storeElement(
398 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>(
399 statePropagatorData, energyData, freeEnergyPerturbationData,
400 globalCommunicationHelper->simulationSignals(),
401 globalCommunicationHelper->nstglobalcomm(), simulator->fplog, simulator->mdlog,
402 simulator->cr, simulator->inputrec, simulator->mdAtoms, simulator->nrnb,
403 simulator->wcycle, simulator->fr, simulator->top_global, simulator->constr));
405 // TODO: Remove this when DD can reduce bonded interactions independently (#3421)
406 auto* castedElement =
407 static_cast<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>*>(
408 vvComputeGlobalsElement);
409 globalCommunicationHelper->setCheckBondedInteractionsCallback(
410 castedElement->getCheckNumberOfBondedInteractionsCallback());
411 builderHelper->storeValue(key, vvComputeGlobalsElement);
412 return vvComputeGlobalsElement;