2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019, 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/nrnb.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/mdlib/md_support.h"
50 #include "gromacs/mdlib/mdatoms.h"
51 #include "gromacs/mdlib/stat.h"
52 #include "gromacs/mdtypes/group.h"
53 #include "gromacs/mdtypes/inputrec.h"
54 #include "gromacs/mdtypes/md_enums.h"
55 #include "gromacs/topology/topology.h"
57 #include "freeenergyperturbationelement.h"
61 template<ComputeGlobalsAlgorithm algorithm>
62 ComputeGlobalsElement<algorithm>::ComputeGlobalsElement(StatePropagatorData* statePropagatorData,
63 EnergyElement* energyElement,
64 FreeEnergyPerturbationElement* freeEnergyPerturbationElement,
65 SimulationSignals* signals,
68 const MDLogger& mdlog,
71 const MDAtoms* mdAtoms,
73 gmx_wallcycle* wcycle,
75 const gmx_mtop_t* global_top,
77 bool hasReadEkinState) :
78 energyReductionStep_(-1),
79 virialReductionStep_(-1),
80 doStopCM_(inputrec->comm_mode != ecmNO),
81 nstcomm_(inputrec->nstcomm),
82 nstglobalcomm_(nstglobalcomm),
83 lastStep_(inputrec->nsteps + inputrec->init_step),
84 initStep_(inputrec->init_step),
85 nullSignaller_(std::make_unique<SimulationSignaller>(nullptr, nullptr, nullptr, false, false)),
86 hasReadEkinState_(hasReadEkinState),
87 totalNumberOfBondedInteractions_(0),
88 shouldCheckNumberOfBondedInteractions_(false),
89 statePropagatorData_(statePropagatorData),
90 energyElement_(energyElement),
91 localTopology_(nullptr),
92 freeEnergyPerturbationElement_(freeEnergyPerturbationElement),
93 vcm_(global_top->groups, *inputrec),
99 top_global_(global_top),
106 reportComRemovalInfo(fplog, vcm_);
107 gstat_ = global_stat_init(inputrec_);
110 template<ComputeGlobalsAlgorithm algorithm>
111 ComputeGlobalsElement<algorithm>::~ComputeGlobalsElement()
113 global_stat_destroy(gstat_);
116 template<ComputeGlobalsAlgorithm algorithm>
117 void ComputeGlobalsElement<algorithm>::elementSetup()
119 GMX_ASSERT(localTopology_, "Setup called before local topology was set.");
121 // Only do initial communication step for one of the velocity-verlet stages
122 if (algorithm == ComputeGlobalsAlgorithm::LeapFrog
123 || algorithm == ComputeGlobalsAlgorithm::VelocityVerletAtFullTimeStep)
125 unsigned int cglo_flags =
126 (CGLO_TEMPERATURE | CGLO_GSTAT
127 | (shouldCheckNumberOfBondedInteractions_ ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
128 | (hasReadEkinState_ ? CGLO_READEKIN : 0));
130 if (algorithm == ComputeGlobalsAlgorithm::VelocityVerletAtFullTimeStep)
132 cglo_flags |= CGLO_PRESSURE | CGLO_CONSTRAINT;
135 const bool stopCM = doStopCM_ && !inputrec_->bContinuation;
137 // To minimize communication, compute_globals computes the COM velocity
138 // and the kinetic energy for the velocities without COM motion removed.
139 // Thus to get the kinetic energy without the COM contribution, we need
140 // to call compute_globals twice.
141 for (int cgloIteration = 0; cgloIteration < (stopCM ? 2 : 1); cgloIteration++)
143 unsigned int cglo_flags_iteration = cglo_flags;
144 if (doStopCM_ && cgloIteration == 0)
146 cglo_flags_iteration |= CGLO_STOPCM;
147 cglo_flags_iteration &= ~CGLO_TEMPERATURE;
150 compute(-1, cglo_flags_iteration, nullSignaller_.get(), false, true);
152 if (cglo_flags_iteration & CGLO_STOPCM)
154 auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
155 // At initialization, do not pass x with acceleration-correction mode
156 // to avoid (incorrect) correction of the initial coordinates.
157 rvec* xPtr = nullptr;
158 if (vcm_.mode != ecmLINEAR_ACCELERATION_CORRECTION)
160 xPtr = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
162 process_and_stopcm_grp(fplog_, &vcm_, *mdAtoms_->mdatoms(), xPtr, v);
163 inc_nrnb(nrnb_, eNR_STOPCM, mdAtoms_->mdatoms()->homenr);
167 // Calculate the initial half step temperature, and save the ekinh_old
168 for (int i = 0; (i < inputrec_->opts.ngtc); i++)
170 copy_mat(energyElement_->ekindata()->tcstat[i].ekinh,
171 energyElement_->ekindata()->tcstat[i].ekinh_old);
176 template<ComputeGlobalsAlgorithm algorithm>
177 void ComputeGlobalsElement<algorithm>::scheduleTask(Step step,
178 Time gmx_unused time,
179 const RegisterRunFunctionPtr& registerRunFunction)
181 const bool needComReduction = doStopCM_ && do_per_step(step, nstcomm_);
182 const bool needGlobalReduction = step == energyReductionStep_ || step == virialReductionStep_
183 || needComReduction || do_per_step(step, nstglobalcomm_);
185 // TODO: CGLO_GSTAT is only used for needToSumEkinhOld_, i.e. to signal that we do or do not
186 // sum the previous kinetic energy. We should simplify / clarify this.
188 if (algorithm == ComputeGlobalsAlgorithm::LeapFrog)
190 // With Leap-Frog we can skip compute_globals at
191 // non-communication steps, but we need to calculate
192 // the kinetic energy one step before communication.
194 // With leap-frog we also need to compute the half-step
195 // kinetic energy at the step before we need to write
196 // the full-step kinetic energy
197 const bool needEkinAtNextStep = (do_per_step(step + 1, nstglobalcomm_) || step + 1 == lastStep_);
199 if (!needGlobalReduction && !needEkinAtNextStep)
204 const bool doEnergy = step == energyReductionStep_;
206 (needGlobalReduction ? CGLO_GSTAT : 0) | (doEnergy ? CGLO_ENERGY : 0)
207 | (needComReduction ? CGLO_STOPCM : 0) | CGLO_TEMPERATURE | CGLO_PRESSURE | CGLO_CONSTRAINT
208 | (shouldCheckNumberOfBondedInteractions_ ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0);
210 // Since we're already communicating at this step, we
211 // can propagate intra-simulation signals. Note that
212 // check_nstglobalcomm has the responsibility for
213 // choosing the value of nstglobalcomm which satisfies
214 // the need of the different signallers.
215 const bool doIntraSimSignal = true;
216 // Disable functionality
217 const bool doInterSimSignal = false;
219 // Make signaller to signal stop / reset / checkpointing signals
220 auto signaller = std::make_shared<SimulationSignaller>(signals_, cr_, nullptr,
221 doInterSimSignal, doIntraSimSignal);
223 (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
224 [this, step, flags, signaller = std::move(signaller)]() {
225 compute(step, flags, signaller.get(), true);
228 else if (algorithm == ComputeGlobalsAlgorithm::VelocityVerletAtFullTimeStep)
230 // For vv, the state at the beginning of the step is positions at time t, velocities at time t - dt/2
231 // The first velocity propagation (+dt/2) therefore actually corresponds to the previous step.
232 // So we need information from the last step in the first half of the integration
233 if (!needGlobalReduction && !do_per_step(step - 1, nstglobalcomm_))
238 const bool doTemperature = step != initStep_ || inputrec_->bContinuation;
239 const bool doEnergy = step == energyReductionStep_;
241 int flags = (needGlobalReduction ? CGLO_GSTAT : 0) | (doEnergy ? CGLO_ENERGY : 0)
242 | (doTemperature ? CGLO_TEMPERATURE : 0) | CGLO_PRESSURE | CGLO_CONSTRAINT
243 | (needComReduction ? CGLO_STOPCM : 0)
244 | (shouldCheckNumberOfBondedInteractions_ ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
247 (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
248 [this, step, flags]() { compute(step, flags, nullSignaller_.get(), false); }));
250 else if (algorithm == ComputeGlobalsAlgorithm::VelocityVerletAfterCoordinateUpdate)
252 // second call to compute_globals for this step
253 if (!needGlobalReduction)
257 int flags = CGLO_GSTAT | CGLO_CONSTRAINT
258 | (shouldCheckNumberOfBondedInteractions_ ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
261 // Since we're already communicating at this step, we
262 // can propagate intra-simulation signals. Note that
263 // check_nstglobalcomm has the responsibility for
264 // choosing the value of nstglobalcomm which satisfies
265 // the need of the different signallers.
266 const bool doIntraSimSignal = true;
267 // Disable functionality
268 const bool doInterSimSignal = false;
270 auto signaller = std::make_shared<SimulationSignaller>(signals_, cr_, nullptr,
271 doInterSimSignal, doIntraSimSignal);
273 (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
274 [this, step, flags, signaller = std::move(signaller)]() {
275 compute(step, flags, signaller.get(), true);
280 template<ComputeGlobalsAlgorithm algorithm>
281 void ComputeGlobalsElement<algorithm>::compute(gmx::Step step,
283 SimulationSignaller* signaller,
287 auto x = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
288 auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
289 auto box = statePropagatorData_->constBox();
290 auto lastbox = useLastBox ? statePropagatorData_->constPreviousBox()
291 : statePropagatorData_->constBox();
293 const real vdwLambda = freeEnergyPerturbationElement_
294 ? freeEnergyPerturbationElement_->constLambdaView()[efptVDW]
297 compute_globals(gstat_, cr_, inputrec_, fr_, energyElement_->ekindata(), x, v, box, vdwLambda,
298 mdAtoms_->mdatoms(), nrnb_, &vcm_, step != -1 ? wcycle_ : nullptr,
299 energyElement_->enerdata(), energyElement_->forceVirial(step),
300 energyElement_->constraintVirial(step), energyElement_->totalVirial(step),
301 energyElement_->pressure(step), energyElement_->muTot(), constr_, signaller, lastbox,
302 &totalNumberOfBondedInteractions_, energyElement_->needToSumEkinhOld(), flags);
303 checkNumberOfBondedInteractions(mdlog_, cr_, totalNumberOfBondedInteractions_, top_global_,
304 localTopology_, x, box, &shouldCheckNumberOfBondedInteractions_);
305 if (flags & CGLO_STOPCM && !isInit)
307 process_and_stopcm_grp(fplog_, &vcm_, *mdAtoms_->mdatoms(), x, v);
308 inc_nrnb(nrnb_, eNR_STOPCM, mdAtoms_->mdatoms()->homenr);
312 template<ComputeGlobalsAlgorithm algorithm>
313 CheckBondedInteractionsCallbackPtr ComputeGlobalsElement<algorithm>::getCheckNumberOfBondedInteractionsCallback()
315 return std::make_unique<CheckBondedInteractionsCallback>(
316 [this]() { needToCheckNumberOfBondedInteractions(); });
319 template<ComputeGlobalsAlgorithm algorithm>
320 void ComputeGlobalsElement<algorithm>::needToCheckNumberOfBondedInteractions()
322 shouldCheckNumberOfBondedInteractions_ = true;
325 template<ComputeGlobalsAlgorithm algorithm>
326 void ComputeGlobalsElement<algorithm>::setTopology(const gmx_localtop_t* top)
328 localTopology_ = top;
331 template<ComputeGlobalsAlgorithm algorithm>
332 SignallerCallbackPtr ComputeGlobalsElement<algorithm>::registerEnergyCallback(EnergySignallerEvent event)
334 if (event == EnergySignallerEvent::EnergyCalculationStep)
336 return std::make_unique<SignallerCallback>(
337 [this](Step step, Time /*unused*/) { energyReductionStep_ = step; });
339 if (event == EnergySignallerEvent::VirialCalculationStep)
341 return std::make_unique<SignallerCallback>(
342 [this](Step step, Time /*unused*/) { virialReductionStep_ = step; });
347 template<ComputeGlobalsAlgorithm algorithm>
348 SignallerCallbackPtr ComputeGlobalsElement<algorithm>::registerTrajectorySignallerCallback(TrajectoryEvent event)
350 if (event == TrajectoryEvent::EnergyWritingStep)
352 return std::make_unique<SignallerCallback>(
353 [this](Step step, Time /*unused*/) { energyReductionStep_ = step; });
358 //! Explicit template instantiation
360 template class ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>;
361 template class ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerletAtFullTimeStep>;
362 template class ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerletAfterCoordinateUpdate>;