src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   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 a velocity-scaling temperature coupling element for
  37  * the modular simulator
  38  *
  39  * \author Pascal Merz <pascal.merz@me.com>
  40  * \ingroup module_modularsimulator
  41  */
  42
  43 #include "gmxpre.h"
  44
  45 #include "velocityscalingtemperaturecoupling.h"
  46
  47 #include <numeric>
  48
  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"
  59
  60 #include "modularsimulator.h"
  61 #include "simulatoralgorithm.h"
  62
  63 namespace gmx
  64 {
  65
  66 /*! \internal
  67  * \brief Data used by the concrete temperature coupling implementations
  68  */
  69 struct TemperatureCouplingData
  70 {
  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;
  81 };
  82
  83 /*! \internal
  84  * \brief Interface for temperature coupling implementations
  85  */
  86 class ITemperatureCouplingImpl
  87 {
  88 public:
  89     //! Allow access to the scaling vectors
  90     virtual void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
  91                                        int numTemperatureGroups) = 0;
  92
  93     /*! \brief Make a temperature control step
  94      *
  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
 100      *
 101      * \return  The temperature coupling integral for the current temperature group
 102      */
 103     [[nodiscard]] virtual real apply(Step                           step,
 104                                      int                            temperatureGroup,
 105                                      real                           currentKineticEnergy,
 106                                      real                           currentTemperature,
 107                                      const TemperatureCouplingData& temperatureCouplingData) = 0;
 108
 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;
 114
 115     //! Standard virtual destructor
 116     virtual ~ITemperatureCouplingImpl() = default;
 117 };
 118
 119 /*! \internal
 120  * \brief Implements v-rescale temperature coupling
 121  */
 122 class VRescaleTemperatureCoupling final : public ITemperatureCouplingImpl
 123 {
 124 public:
 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
 131     {
 132         if (!(temperatureCouplingData.couplingTime[temperatureGroup] >= 0
 133               && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0
 134               && currentKineticEnergy > 0))
 135         {
 136             lambdaStartVelocities_[temperatureGroup] = 1.0;
 137             return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
 138         }
 139
 140         const real referenceKineticEnergy =
 141                 0.5 * temperatureCouplingData.referenceTemperature[temperatureGroup] * BOLTZ
 142                 * temperatureCouplingData.numDegreesOfFreedom[temperatureGroup];
 143
 144         const real newKineticEnergy =
 145                 vrescale_resamplekin(currentKineticEnergy,
 146                                      referenceKineticEnergy,
 147                                      temperatureCouplingData.numDegreesOfFreedom[temperatureGroup],
 148                                      temperatureCouplingData.couplingTime[temperatureGroup]
 149                                              / temperatureCouplingData.couplingTimeStep,
 150                                      step,
 151                                      seed_);
 152
 153         // Analytically newKineticEnergy >= 0, but we check for rounding errors
 154         if (newKineticEnergy <= 0)
 155         {
 156             lambdaStartVelocities_[temperatureGroup] = 0.0;
 157         }
 158         else
 159         {
 160             lambdaStartVelocities_[temperatureGroup] = std::sqrt(newKineticEnergy / currentKineticEnergy);
 161         }
 162
 163         if (debug)
 164         {
 165             fprintf(debug,
 166                     "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
 167                     temperatureGroup,
 168                     referenceKineticEnergy,
 169                     currentKineticEnergy,
 170                     newKineticEnergy,
 171                     lambdaStartVelocities_[temperatureGroup]);
 172         }
 173
 174         return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup]
 175                - (newKineticEnergy - currentKineticEnergy);
 176     }
 177
 178     //! Connect with propagator - v-rescale only scales start step velocities
 179     void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
 180                                int                                   numTemperatureGroups) override
 181     {
 182         connectionData.setNumVelocityScalingVariables(numTemperatureGroups);
 183         lambdaStartVelocities_ = connectionData.getViewOnVelocityScaling();
 184     }
 185
 186     //! No data to write to checkpoint
 187     void writeCheckpoint(std::optional<WriteCheckpointData> gmx_unused checkpointData,
 188                          const t_commrec gmx_unused* cr) override
 189     {
 190     }
 191     //! No data to read from checkpoints
 192     void readCheckpoint(std::optional<ReadCheckpointData> gmx_unused checkpointData,
 193                         const t_commrec gmx_unused* cr) override
 194     {
 195     }
 196
 197     //! Constructor
 198     VRescaleTemperatureCoupling(int64_t seed) : seed_(seed) {}
 199
 200 private:
 201     //! The random seed
 202     const int64_t seed_;
 203
 204     //! View on the scaling factor of the propagator (pre-step velocities)
 205     ArrayRef<real> lambdaStartVelocities_;
 206 };
 207
 208 /*! \internal
 209  * \brief Implements Berendsen temperature coupling
 210  */
 211 class BerendsenTemperatureCoupling final : public ITemperatureCouplingImpl
 212 {
 213 public:
 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
 220     {
 221         if (!(temperatureCouplingData.couplingTime[temperatureGroup] >= 0
 222               && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0
 223               && currentKineticEnergy > 0))
 224         {
 225             lambdaStartVelocities_[temperatureGroup] = 1.0;
 226             return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
 227         }
 228
 229         real lambda =
 230                 std::sqrt(1.0
 231                           + (temperatureCouplingData.couplingTimeStep
 232                              / temperatureCouplingData.couplingTime[temperatureGroup])
 233                                     * (temperatureCouplingData.referenceTemperature[temperatureGroup] / currentTemperature
 234                                        - 1.0));
 235         lambdaStartVelocities_[temperatureGroup] =
 236                 std::max<real>(std::min<real>(lambda, 1.25_real), 0.8_real);
 237         if (debug)
 238         {
 239             fprintf(debug,
 240                     "TC: group %d: T: %g, Lambda: %g\n",
 241                     temperatureGroup,
 242                     currentTemperature,
 243                     lambdaStartVelocities_[temperatureGroup]);
 244         }
 245         return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup]
 246                - (lambdaStartVelocities_[temperatureGroup] * lambdaStartVelocities_[temperatureGroup]
 247                   - 1) * currentKineticEnergy;
 248     }
 249
 250     //! Connect with propagator - Berendsen only scales start step velocities
 251     void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
 252                                int                                   numTemperatureGroups) override
 253     {
 254         connectionData.setNumVelocityScalingVariables(numTemperatureGroups);
 255         lambdaStartVelocities_ = connectionData.getViewOnVelocityScaling();
 256     }
 257
 258     //! No data to write to checkpoint
 259     void writeCheckpoint(std::optional<WriteCheckpointData> gmx_unused checkpointData,
 260                          const t_commrec gmx_unused* cr) override
 261     {
 262     }
 263     //! No data to read from checkpoints
 264     void readCheckpoint(std::optional<ReadCheckpointData> gmx_unused checkpointData,
 265                         const t_commrec gmx_unused* cr) override
 266     {
 267     }
 268
 269 private:
 270     //! View on the scaling factor of the propagator (pre-step velocities)
 271     ArrayRef<real> lambdaStartVelocities_;
 272 };
 273
 274 VelocityScalingTemperatureCoupling::VelocityScalingTemperatureCoupling(
 275         int                               nstcouple,
 276         int                               offset,
 277         UseFullStepKE                     useFullStepKE,
 278         ReportPreviousStepConservedEnergy reportPreviousConservedEnergy,
 279         int64_t                           seed,
 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),
 288     offset_(offset),
 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)
 298 {
 299     if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
 300     {
 301         temperatureCouplingIntegralPreviousStep_ = temperatureCouplingIntegral_;
 302     }
 303     energyData->setVelocityScalingTemperatureCoupling(this);
 304     if (couplingType == TemperatureCoupling::VRescale)
 305     {
 306         temperatureCouplingImpl_ = std::make_unique<VRescaleTemperatureCoupling>(seed);
 307     }
 308     else if (couplingType == TemperatureCoupling::Berendsen)
 309     {
 310         temperatureCouplingImpl_ = std::make_unique<BerendsenTemperatureCoupling>();
 311     }
 312     else
 313     {
 314         throw NotImplementedError("Temperature coupling " + std::string(enumValueToString(couplingType))
 315                                   + " is not implemented for modular simulator.");
 316     }
 317 }
 318
 319 void VelocityScalingTemperatureCoupling::connectWithPropagator(const PropagatorThermostatConnection& connectionData)
 320 {
 321     temperatureCouplingImpl_->connectWithPropagator(connectionData, numTemperatureGroups_);
 322     propagatorCallback_ = connectionData.getVelocityScalingCallback();
 323 }
 324
 325 void VelocityScalingTemperatureCoupling::elementSetup()
 326 {
 327     if (!propagatorCallback_)
 328     {
 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 "
 333                 "object.");
 334     }
 335 }
 336
 337 void VelocityScalingTemperatureCoupling::scheduleTask(Step step,
 338                                                       Time gmx_unused            time,
 339                                                       const RegisterRunFunction& registerRunFunction)
 340 {
 341     /* The thermostat will need a valid kinetic energy when it is running.
 342      * Currently, computeGlobalCommunicationPeriod() is making sure this
 343      * happens on time.
 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.
 347      *
 348      */
 349     if (do_per_step(step + nstcouple_ + offset_, nstcouple_))
 350     {
 351         // do T-coupling this step
 352         registerRunFunction([this, step]() { setLambda(step); });
 353
 354         // Let propagator know that we want to do T-coupling
 355         propagatorCallback_(step);
 356     }
 357 }
 358
 359 void VelocityScalingTemperatureCoupling::setLambda(Step step)
 360 {
 361     // if we report the previous energy, calculate before the step
 362     if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
 363     {
 364         temperatureCouplingIntegralPreviousStep_ = temperatureCouplingIntegral_;
 365     }
 366
 367     const auto*             ekind          = energyData_->ekindata();
 368     TemperatureCouplingData thermostatData = {
 369         couplingTimeStep_, referenceTemperature_, couplingTime_, numDegreesOfFreedom_, temperatureCouplingIntegral_
 370     };
 371
 372     for (int temperatureGroup = 0; (temperatureGroup < numTemperatureGroups_); temperatureGroup++)
 373     {
 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;
 380
 381         temperatureCouplingIntegral_[temperatureGroup] = temperatureCouplingImpl_->apply(
 382                 step, temperatureGroup, currentKineticEnergy, currentTemperature, thermostatData);
 383     }
 384 }
 385
 386 namespace
 387 {
 388 /*!
 389  * \brief Enum describing the contents VelocityScalingTemperatureCoupling writes to modular checkpoint
 390  *
 391  * When changing the checkpoint content, add a new element just above Count, and adjust the
 392  * checkpoint functionality.
 393  */
 394 enum class CheckpointVersion
 395 {
 396     Base, //!< First version of modular checkpointing
 397     Count //!< Number of entries. Add new versions right above this!
 398 };
 399 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
 400 } // namespace
 401
 402 template<CheckpointDataOperation operation>
 403 void VelocityScalingTemperatureCoupling::doCheckpointData(CheckpointData<operation>* checkpointData)
 404 {
 405     checkpointVersion(checkpointData, "VRescaleThermostat version", c_currentVersion);
 406
 407     checkpointData->arrayRef("thermostat integral",
 408                              makeCheckpointArrayRef<operation>(temperatureCouplingIntegral_));
 409 }
 410
 411 void VelocityScalingTemperatureCoupling::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
 412                                                              const t_commrec*                   cr)
 413 {
 414     if (MASTER(cr))
 415     {
 416         doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
 417     }
 418     temperatureCouplingImpl_->writeCheckpoint(
 419             checkpointData
 420                     ? std::make_optional(checkpointData->subCheckpointData("thermostat impl"))
 421                     : std::nullopt,
 422             cr);
 423 }
 424
 425 void VelocityScalingTemperatureCoupling::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
 426                                                                 const t_commrec* cr)
 427 {
 428     if (MASTER(cr))
 429     {
 430         doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
 431     }
 432     if (DOMAINDECOMP(cr))
 433     {
 434         dd_bcast(cr->dd,
 435                  ssize(temperatureCouplingIntegral_) * int(sizeof(double)),
 436                  temperatureCouplingIntegral_.data());
 437     }
 438     temperatureCouplingImpl_->readCheckpoint(
 439             checkpointData
 440                     ? std::make_optional(checkpointData->subCheckpointData("thermostat impl"))
 441                     : std::nullopt,
 442             cr);
 443 }
 444
 445 const std::string& VelocityScalingTemperatureCoupling::clientID()
 446 {
 447     return identifier_;
 448 }
 449
 450 real VelocityScalingTemperatureCoupling::conservedEnergyContribution() const
 451 {
 452     if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
 453     {
 454         return std::accumulate(temperatureCouplingIntegralPreviousStep_.begin(),
 455                                temperatureCouplingIntegralPreviousStep_.end(),
 456                                0.0);
 457     }
 458     else
 459     {
 460         return std::accumulate(
 461                 temperatureCouplingIntegral_.begin(), temperatureCouplingIntegral_.end(), 0.0);
 462     }
 463 }
 464
 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,
 472         int                                   offset,
 473         UseFullStepKE                         useFullStepKE,
 474         ReportPreviousStepConservedEnergy     reportPreviousStepConservedEnergy)
 475 {
 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,
 479             offset,
 480             useFullStepKE,
 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,
 488             energyData,
 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);
 494     });
 495     return element;
 496 }
 497
 498 } // namespace gmx