2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020,2021, 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 a velocity-scaling temperature coupling element for
37 * the modular simulator
39 * \author Pascal Merz <pascal.merz@me.com>
40 * \ingroup module_modularsimulator
45 #include "velocityscalingtemperaturecoupling.h"
49 #include "gromacs/domdec/domdec_network.h"
50 #include "gromacs/math/units.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdlib/coupling.h"
53 #include "gromacs/mdlib/stat.h"
54 #include "gromacs/mdtypes/checkpointdata.h"
55 #include "gromacs/mdtypes/commrec.h"
56 #include "gromacs/mdtypes/group.h"
57 #include "gromacs/mdtypes/inputrec.h"
58 #include "gromacs/utility/fatalerror.h"
60 #include "modularsimulator.h"
61 #include "simulatoralgorithm.h"
67 * \brief Data used by the concrete temperature coupling implementations
69 struct TemperatureCouplingData
71 //! The coupling time step - simulation time step x nstcouple_
72 const double couplingTimeStep;
73 //! Coupling temperature per group
74 ArrayRef<const real> referenceTemperature;
75 //! Coupling time per group
76 ArrayRef<const real> couplingTime;
77 //! Number of degrees of freedom per group
78 ArrayRef<const real> numDegreesOfFreedom;
79 //! Work exerted by thermostat per group
80 ArrayRef<const double> temperatureCouplingIntegral;
84 * \brief Interface for temperature coupling implementations
86 class ITemperatureCouplingImpl
89 //! Allow access to the scaling vectors
90 virtual void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
91 int numTemperatureGroups) = 0;
93 /*! \brief Make a temperature control step
95 * \param step The current step
96 * \param temperatureGroup The current temperature group
97 * \param currentKineticEnergy The kinetic energy of the temperature group
98 * \param currentTemperature The temperature of the temperature group
99 * \param temperatureCouplingData Access to general temperature coupling data
101 * \return The temperature coupling integral for the current temperature group
103 [[nodiscard]] virtual real apply(Step step,
104 int temperatureGroup,
105 real currentKineticEnergy,
106 real currentTemperature,
107 const TemperatureCouplingData& temperatureCouplingData) = 0;
109 //! Write private data to checkpoint
110 virtual void writeCheckpoint(std::optional<WriteCheckpointData> checkpointData,
111 const t_commrec* cr) = 0;
112 //! Read private data from checkpoint
113 virtual void readCheckpoint(std::optional<ReadCheckpointData> checkpointData, const t_commrec* cr) = 0;
115 //! Standard virtual destructor
116 virtual ~ITemperatureCouplingImpl() = default;
120 * \brief Implements v-rescale temperature coupling
122 class VRescaleTemperatureCoupling final : public ITemperatureCouplingImpl
125 //! Apply the v-rescale temperature control
126 real apply(Step step,
127 int temperatureGroup,
128 real currentKineticEnergy,
129 real gmx_unused currentTemperature,
130 const TemperatureCouplingData& temperatureCouplingData) override
132 if (!(temperatureCouplingData.couplingTime[temperatureGroup] >= 0
133 && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0
134 && currentKineticEnergy > 0))
136 lambdaStartVelocities_[temperatureGroup] = 1.0;
137 return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
140 const real referenceKineticEnergy =
141 0.5 * temperatureCouplingData.referenceTemperature[temperatureGroup] * gmx::c_boltz
142 * temperatureCouplingData.numDegreesOfFreedom[temperatureGroup];
144 const real newKineticEnergy =
145 vrescale_resamplekin(currentKineticEnergy,
146 referenceKineticEnergy,
147 temperatureCouplingData.numDegreesOfFreedom[temperatureGroup],
148 temperatureCouplingData.couplingTime[temperatureGroup]
149 / temperatureCouplingData.couplingTimeStep,
153 // Analytically newKineticEnergy >= 0, but we check for rounding errors
154 if (newKineticEnergy <= 0)
156 lambdaStartVelocities_[temperatureGroup] = 0.0;
160 lambdaStartVelocities_[temperatureGroup] = std::sqrt(newKineticEnergy / currentKineticEnergy);
166 "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
168 referenceKineticEnergy,
169 currentKineticEnergy,
171 lambdaStartVelocities_[temperatureGroup]);
174 return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup]
175 - (newKineticEnergy - currentKineticEnergy);
178 //! Connect with propagator - v-rescale only scales start step velocities
179 void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
180 int numTemperatureGroups) override
182 connectionData.setNumVelocityScalingVariables(numTemperatureGroups);
183 lambdaStartVelocities_ = connectionData.getViewOnVelocityScaling();
186 //! No data to write to checkpoint
187 void writeCheckpoint(std::optional<WriteCheckpointData> gmx_unused checkpointData,
188 const t_commrec gmx_unused* cr) override
191 //! No data to read from checkpoints
192 void readCheckpoint(std::optional<ReadCheckpointData> gmx_unused checkpointData,
193 const t_commrec gmx_unused* cr) override
198 VRescaleTemperatureCoupling(int64_t seed) : seed_(seed) {}
204 //! View on the scaling factor of the propagator (pre-step velocities)
205 ArrayRef<real> lambdaStartVelocities_;
209 * \brief Implements Berendsen temperature coupling
211 class BerendsenTemperatureCoupling final : public ITemperatureCouplingImpl
214 //! Apply the v-rescale temperature control
215 real apply(Step gmx_unused step,
216 int temperatureGroup,
217 real currentKineticEnergy,
218 real currentTemperature,
219 const TemperatureCouplingData& temperatureCouplingData) override
221 if (!(temperatureCouplingData.couplingTime[temperatureGroup] >= 0
222 && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0
223 && currentKineticEnergy > 0))
225 lambdaStartVelocities_[temperatureGroup] = 1.0;
226 return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
231 + (temperatureCouplingData.couplingTimeStep
232 / temperatureCouplingData.couplingTime[temperatureGroup])
233 * (temperatureCouplingData.referenceTemperature[temperatureGroup] / currentTemperature
235 lambdaStartVelocities_[temperatureGroup] =
236 std::max<real>(std::min<real>(lambda, 1.25_real), 0.8_real);
240 "TC: group %d: T: %g, Lambda: %g\n",
243 lambdaStartVelocities_[temperatureGroup]);
245 return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup]
246 - (lambdaStartVelocities_[temperatureGroup] * lambdaStartVelocities_[temperatureGroup]
247 - 1) * currentKineticEnergy;
250 //! Connect with propagator - Berendsen only scales start step velocities
251 void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
252 int numTemperatureGroups) override
254 connectionData.setNumVelocityScalingVariables(numTemperatureGroups);
255 lambdaStartVelocities_ = connectionData.getViewOnVelocityScaling();
258 //! No data to write to checkpoint
259 void writeCheckpoint(std::optional<WriteCheckpointData> gmx_unused checkpointData,
260 const t_commrec gmx_unused* cr) override
263 //! No data to read from checkpoints
264 void readCheckpoint(std::optional<ReadCheckpointData> gmx_unused checkpointData,
265 const t_commrec gmx_unused* cr) override
270 //! View on the scaling factor of the propagator (pre-step velocities)
271 ArrayRef<real> lambdaStartVelocities_;
274 VelocityScalingTemperatureCoupling::VelocityScalingTemperatureCoupling(
277 UseFullStepKE useFullStepKE,
278 ReportPreviousStepConservedEnergy reportPreviousConservedEnergy,
280 int numTemperatureGroups,
281 double couplingTimeStep,
282 const real* referenceTemperature,
283 const real* couplingTime,
284 const real* numDegreesOfFreedom,
285 EnergyData* energyData,
286 TemperatureCoupling couplingType) :
287 nstcouple_(nstcouple),
289 useFullStepKE_(useFullStepKE),
290 reportPreviousConservedEnergy_(reportPreviousConservedEnergy),
291 numTemperatureGroups_(numTemperatureGroups),
292 couplingTimeStep_(couplingTimeStep),
293 referenceTemperature_(referenceTemperature, referenceTemperature + numTemperatureGroups),
294 couplingTime_(couplingTime, couplingTime + numTemperatureGroups),
295 numDegreesOfFreedom_(numDegreesOfFreedom, numDegreesOfFreedom + numTemperatureGroups),
296 temperatureCouplingIntegral_(numTemperatureGroups, 0.0),
297 energyData_(energyData)
299 if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
301 temperatureCouplingIntegralPreviousStep_ = temperatureCouplingIntegral_;
303 energyData->setVelocityScalingTemperatureCoupling(this);
304 if (couplingType == TemperatureCoupling::VRescale)
306 temperatureCouplingImpl_ = std::make_unique<VRescaleTemperatureCoupling>(seed);
308 else if (couplingType == TemperatureCoupling::Berendsen)
310 temperatureCouplingImpl_ = std::make_unique<BerendsenTemperatureCoupling>();
314 throw NotImplementedError("Temperature coupling " + std::string(enumValueToString(couplingType))
315 + " is not implemented for modular simulator.");
319 void VelocityScalingTemperatureCoupling::connectWithPropagator(const PropagatorThermostatConnection& connectionData)
321 temperatureCouplingImpl_->connectWithPropagator(connectionData, numTemperatureGroups_);
322 propagatorCallback_ = connectionData.getVelocityScalingCallback();
325 void VelocityScalingTemperatureCoupling::elementSetup()
327 if (!propagatorCallback_)
329 throw MissingElementConnectionError(
330 "Velocity scaling temperature coupling was not connected to a propagator.\n"
331 "Connection to a propagator element is needed to scale the velocities.\n"
332 "Use connectWithPropagator(...) before building the ModularSimulatorAlgorithm "
337 void VelocityScalingTemperatureCoupling::scheduleTask(Step step,
338 Time gmx_unused time,
339 const RegisterRunFunction& registerRunFunction)
341 /* The thermostat will need a valid kinetic energy when it is running.
342 * Currently, computeGlobalCommunicationPeriod() is making sure this
344 * TODO: Once we're switching to a new global communication scheme, we
345 * will want the thermostat to signal that global reduction
346 * of the kinetic energy is needed.
349 if (do_per_step(step + nstcouple_ + offset_, nstcouple_))
351 // do T-coupling this step
352 registerRunFunction([this, step]() { setLambda(step); });
354 // Let propagator know that we want to do T-coupling
355 propagatorCallback_(step);
359 void VelocityScalingTemperatureCoupling::setLambda(Step step)
361 // if we report the previous energy, calculate before the step
362 if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
364 temperatureCouplingIntegralPreviousStep_ = temperatureCouplingIntegral_;
367 const auto* ekind = energyData_->ekindata();
368 TemperatureCouplingData thermostatData = {
369 couplingTimeStep_, referenceTemperature_, couplingTime_, numDegreesOfFreedom_, temperatureCouplingIntegral_
372 for (int temperatureGroup = 0; (temperatureGroup < numTemperatureGroups_); temperatureGroup++)
374 const real currentKineticEnergy = useFullStepKE_ == UseFullStepKE::Yes
375 ? trace(ekind->tcstat[temperatureGroup].ekinf)
376 : trace(ekind->tcstat[temperatureGroup].ekinh);
377 const real currentTemperature = useFullStepKE_ == UseFullStepKE::Yes
378 ? ekind->tcstat[temperatureGroup].T
379 : ekind->tcstat[temperatureGroup].Th;
381 temperatureCouplingIntegral_[temperatureGroup] = temperatureCouplingImpl_->apply(
382 step, temperatureGroup, currentKineticEnergy, currentTemperature, thermostatData);
389 * \brief Enum describing the contents VelocityScalingTemperatureCoupling writes to modular checkpoint
391 * When changing the checkpoint content, add a new element just above Count, and adjust the
392 * checkpoint functionality.
394 enum class CheckpointVersion
396 Base, //!< First version of modular checkpointing
397 Count //!< Number of entries. Add new versions right above this!
399 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
402 template<CheckpointDataOperation operation>
403 void VelocityScalingTemperatureCoupling::doCheckpointData(CheckpointData<operation>* checkpointData)
405 checkpointVersion(checkpointData, "VRescaleThermostat version", c_currentVersion);
407 checkpointData->arrayRef("thermostat integral",
408 makeCheckpointArrayRef<operation>(temperatureCouplingIntegral_));
411 void VelocityScalingTemperatureCoupling::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
416 doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
418 temperatureCouplingImpl_->writeCheckpoint(
420 ? std::make_optional(checkpointData->subCheckpointData("thermostat impl"))
425 void VelocityScalingTemperatureCoupling::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
430 doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
432 if (DOMAINDECOMP(cr))
435 ssize(temperatureCouplingIntegral_) * int(sizeof(double)),
436 temperatureCouplingIntegral_.data());
438 temperatureCouplingImpl_->readCheckpoint(
440 ? std::make_optional(checkpointData->subCheckpointData("thermostat impl"))
445 const std::string& VelocityScalingTemperatureCoupling::clientID()
450 real VelocityScalingTemperatureCoupling::conservedEnergyContribution() const
452 if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
454 return std::accumulate(temperatureCouplingIntegralPreviousStep_.begin(),
455 temperatureCouplingIntegralPreviousStep_.end(),
460 return std::accumulate(
461 temperatureCouplingIntegral_.begin(), temperatureCouplingIntegral_.end(), 0.0);
465 ISimulatorElement* VelocityScalingTemperatureCoupling::getElementPointerImpl(
466 LegacySimulatorData* legacySimulatorData,
467 ModularSimulatorAlgorithmBuilderHelper* builderHelper,
468 StatePropagatorData gmx_unused* statePropagatorData,
469 EnergyData gmx_unused* energyData,
470 FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
471 GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
473 UseFullStepKE useFullStepKE,
474 ReportPreviousStepConservedEnergy reportPreviousStepConservedEnergy)
476 // Element is now owned by the caller of this method, who will handle lifetime (see ModularSimulatorAlgorithm)
477 auto* element = builderHelper->storeElement(std::make_unique<VelocityScalingTemperatureCoupling>(
478 legacySimulatorData->inputrec->nsttcouple,
481 reportPreviousStepConservedEnergy,
482 legacySimulatorData->inputrec->ld_seed,
483 legacySimulatorData->inputrec->opts.ngtc,
484 legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nsttcouple,
485 legacySimulatorData->inputrec->opts.ref_t,
486 legacySimulatorData->inputrec->opts.tau_t,
487 legacySimulatorData->inputrec->opts.nrdf,
489 legacySimulatorData->inputrec->etc));
490 auto* thermostat = static_cast<VelocityScalingTemperatureCoupling*>(element);
491 // Capturing pointer is safe because lifetime is handled by caller
492 builderHelper->registerThermostat([thermostat](const PropagatorThermostatConnection& connection) {
493 thermostat->connectWithPropagator(connection);