src/gromacs/modularsimulator/velocityscalingtemperaturecoupling.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 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, referenceKineticEnergy,
 146                                      temperatureCouplingData.numDegreesOfFreedom[temperatureGroup],
 147                                      temperatureCouplingData.couplingTime[temperatureGroup]
 148                                              / temperatureCouplingData.couplingTimeStep,
 149                                      step, seed_);
 150
 151         // Analytically newKineticEnergy >= 0, but we check for rounding errors
 152         if (newKineticEnergy <= 0)
 153         {
 154             lambdaStartVelocities_[temperatureGroup] = 0.0;
 155         }
 156         else
 157         {
 158             lambdaStartVelocities_[temperatureGroup] = std::sqrt(newKineticEnergy / currentKineticEnergy);
 159         }
 160
 161         if (debug)
 162         {
 163             fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n", temperatureGroup,
 164                     referenceKineticEnergy, currentKineticEnergy, newKineticEnergy,
 165                     lambdaStartVelocities_[temperatureGroup]);
 166         }
 167
 168         return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup]
 169                - (newKineticEnergy - currentKineticEnergy);
 170     }
 171
 172     //! Connect with propagator - v-rescale only scales start step velocities
 173     void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
 174                                int                                   numTemperatureGroups) override
 175     {
 176         connectionData.setNumVelocityScalingVariables(numTemperatureGroups);
 177         lambdaStartVelocities_ = connectionData.getViewOnVelocityScaling();
 178     }
 179
 180     //! No data to write to checkpoint
 181     void writeCheckpoint(std::optional<WriteCheckpointData> gmx_unused checkpointData,
 182                          const t_commrec gmx_unused* cr) override
 183     {
 184     }
 185     //! No data to read from checkpoints
 186     void readCheckpoint(std::optional<ReadCheckpointData> gmx_unused checkpointData,
 187                         const t_commrec gmx_unused* cr) override
 188     {
 189     }
 190
 191     //! Constructor
 192     VRescaleTemperatureCoupling(int64_t seed) : seed_(seed) {}
 193
 194 private:
 195     //! The random seed
 196     const int64_t seed_;
 197
 198     //! View on the scaling factor of the propagator (pre-step velocities)
 199     ArrayRef<real> lambdaStartVelocities_;
 200 };
 201
 202 /*! \internal
 203  * \brief Implements Berendsen temperature coupling
 204  */
 205 class BerendsenTemperatureCoupling final : public ITemperatureCouplingImpl
 206 {
 207 public:
 208     //! Apply the v-rescale temperature control
 209     real apply(Step gmx_unused                step,
 210                int                            temperatureGroup,
 211                real                           currentKineticEnergy,
 212                real                           currentTemperature,
 213                const TemperatureCouplingData& temperatureCouplingData) override
 214     {
 215         if (!(temperatureCouplingData.couplingTime[temperatureGroup] >= 0
 216               && temperatureCouplingData.numDegreesOfFreedom[temperatureGroup] > 0
 217               && currentKineticEnergy > 0))
 218         {
 219             lambdaStartVelocities_[temperatureGroup] = 1.0;
 220             return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup];
 221         }
 222
 223         real lambda =
 224                 std::sqrt(1.0
 225                           + (temperatureCouplingData.couplingTimeStep
 226                              / temperatureCouplingData.couplingTime[temperatureGroup])
 227                                     * (temperatureCouplingData.referenceTemperature[temperatureGroup] / currentTemperature
 228                                        - 1.0));
 229         lambdaStartVelocities_[temperatureGroup] =
 230                 std::max<real>(std::min<real>(lambda, 1.25_real), 0.8_real);
 231         if (debug)
 232         {
 233             fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", temperatureGroup,
 234                     currentTemperature, lambdaStartVelocities_[temperatureGroup]);
 235         }
 236         return temperatureCouplingData.temperatureCouplingIntegral[temperatureGroup]
 237                - (lambdaStartVelocities_[temperatureGroup] * lambdaStartVelocities_[temperatureGroup]
 238                   - 1) * currentKineticEnergy;
 239     }
 240
 241     //! Connect with propagator - Berendsen only scales start step velocities
 242     void connectWithPropagator(const PropagatorThermostatConnection& connectionData,
 243                                int                                   numTemperatureGroups) override
 244     {
 245         connectionData.setNumVelocityScalingVariables(numTemperatureGroups);
 246         lambdaStartVelocities_ = connectionData.getViewOnVelocityScaling();
 247     }
 248
 249     //! No data to write to checkpoint
 250     void writeCheckpoint(std::optional<WriteCheckpointData> gmx_unused checkpointData,
 251                          const t_commrec gmx_unused* cr) override
 252     {
 253     }
 254     //! No data to read from checkpoints
 255     void readCheckpoint(std::optional<ReadCheckpointData> gmx_unused checkpointData,
 256                         const t_commrec gmx_unused* cr) override
 257     {
 258     }
 259
 260 private:
 261     //! View on the scaling factor of the propagator (pre-step velocities)
 262     ArrayRef<real> lambdaStartVelocities_;
 263 };
 264
 265 VelocityScalingTemperatureCoupling::VelocityScalingTemperatureCoupling(
 266         int                               nstcouple,
 267         int                               offset,
 268         UseFullStepKE                     useFullStepKE,
 269         ReportPreviousStepConservedEnergy reportPreviousConservedEnergy,
 270         int64_t                           seed,
 271         int                               numTemperatureGroups,
 272         double                            couplingTimeStep,
 273         const real*                       referenceTemperature,
 274         const real*                       couplingTime,
 275         const real*                       numDegreesOfFreedom,
 276         EnergyData*                       energyData,
 277         TemperatureCouplingType           couplingType) :
 278     nstcouple_(nstcouple),
 279     offset_(offset),
 280     useFullStepKE_(useFullStepKE),
 281     reportPreviousConservedEnergy_(reportPreviousConservedEnergy),
 282     numTemperatureGroups_(numTemperatureGroups),
 283     couplingTimeStep_(couplingTimeStep),
 284     referenceTemperature_(referenceTemperature, referenceTemperature + numTemperatureGroups),
 285     couplingTime_(couplingTime, couplingTime + numTemperatureGroups),
 286     numDegreesOfFreedom_(numDegreesOfFreedom, numDegreesOfFreedom + numTemperatureGroups),
 287     temperatureCouplingIntegral_(numTemperatureGroups, 0.0),
 288     energyData_(energyData)
 289 {
 290     if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
 291     {
 292         temperatureCouplingIntegralPreviousStep_ = temperatureCouplingIntegral_;
 293     }
 294     energyData->setVelocityScalingTemperatureCoupling(this);
 295     if (couplingType == etcVRESCALE)
 296     {
 297         temperatureCouplingImpl_ = std::make_unique<VRescaleTemperatureCoupling>(seed);
 298     }
 299     else if (couplingType == etcBERENDSEN)
 300     {
 301         temperatureCouplingImpl_ = std::make_unique<BerendsenTemperatureCoupling>();
 302     }
 303     else
 304     {
 305         throw NotImplementedError("Temperature coupling " + std::string(ETCOUPLTYPE(couplingType))
 306                                   + " is not implemented for modular simulator.");
 307     }
 308 }
 309
 310 void VelocityScalingTemperatureCoupling::connectWithPropagator(const PropagatorThermostatConnection& connectionData)
 311 {
 312     temperatureCouplingImpl_->connectWithPropagator(connectionData, numTemperatureGroups_);
 313     propagatorCallback_ = connectionData.getVelocityScalingCallback();
 314 }
 315
 316 void VelocityScalingTemperatureCoupling::elementSetup()
 317 {
 318     if (!propagatorCallback_)
 319     {
 320         throw MissingElementConnectionError(
 321                 "Velocity scaling temperature coupling was not connected to a propagator.\n"
 322                 "Connection to a propagator element is needed to scale the velocities.\n"
 323                 "Use connectWithPropagator(...) before building the ModularSimulatorAlgorithm "
 324                 "object.");
 325     }
 326 }
 327
 328 void VelocityScalingTemperatureCoupling::scheduleTask(Step step,
 329                                                       Time gmx_unused            time,
 330                                                       const RegisterRunFunction& registerRunFunction)
 331 {
 332     /* The thermostat will need a valid kinetic energy when it is running.
 333      * Currently, computeGlobalCommunicationPeriod() is making sure this
 334      * happens on time.
 335      * TODO: Once we're switching to a new global communication scheme, we
 336      *       will want the thermostat to signal that global reduction
 337      *       of the kinetic energy is needed.
 338      *
 339      */
 340     if (do_per_step(step + nstcouple_ + offset_, nstcouple_))
 341     {
 342         // do T-coupling this step
 343         registerRunFunction([this, step]() { setLambda(step); });
 344
 345         // Let propagator know that we want to do T-coupling
 346         propagatorCallback_(step);
 347     }
 348 }
 349
 350 void VelocityScalingTemperatureCoupling::setLambda(Step step)
 351 {
 352     // if we report the previous energy, calculate before the step
 353     if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
 354     {
 355         temperatureCouplingIntegralPreviousStep_ = temperatureCouplingIntegral_;
 356     }
 357
 358     const auto*             ekind          = energyData_->ekindata();
 359     TemperatureCouplingData thermostatData = { couplingTimeStep_, referenceTemperature_, couplingTime_,
 360                                                numDegreesOfFreedom_, temperatureCouplingIntegral_ };
 361
 362     for (int temperatureGroup = 0; (temperatureGroup < numTemperatureGroups_); temperatureGroup++)
 363     {
 364         const real currentKineticEnergy = useFullStepKE_ == UseFullStepKE::Yes
 365                                                   ? trace(ekind->tcstat[temperatureGroup].ekinf)
 366                                                   : trace(ekind->tcstat[temperatureGroup].ekinh);
 367         const real currentTemperature = useFullStepKE_ == UseFullStepKE::Yes
 368                                                 ? ekind->tcstat[temperatureGroup].T
 369                                                 : ekind->tcstat[temperatureGroup].Th;
 370
 371         temperatureCouplingIntegral_[temperatureGroup] = temperatureCouplingImpl_->apply(
 372                 step, temperatureGroup, currentKineticEnergy, currentTemperature, thermostatData);
 373     }
 374 }
 375
 376 namespace
 377 {
 378 /*!
 379  * \brief Enum describing the contents VelocityScalingTemperatureCoupling writes to modular checkpoint
 380  *
 381  * When changing the checkpoint content, add a new element just above Count, and adjust the
 382  * checkpoint functionality.
 383  */
 384 enum class CheckpointVersion
 385 {
 386     Base, //!< First version of modular checkpointing
 387     Count //!< Number of entries. Add new versions right above this!
 388 };
 389 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
 390 } // namespace
 391
 392 template<CheckpointDataOperation operation>
 393 void VelocityScalingTemperatureCoupling::doCheckpointData(CheckpointData<operation>* checkpointData)
 394 {
 395     checkpointVersion(checkpointData, "VRescaleThermostat version", c_currentVersion);
 396
 397     checkpointData->arrayRef("thermostat integral",
 398                              makeCheckpointArrayRef<operation>(temperatureCouplingIntegral_));
 399 }
 400
 401 void VelocityScalingTemperatureCoupling::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
 402                                                              const t_commrec*                   cr)
 403 {
 404     if (MASTER(cr))
 405     {
 406         doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
 407     }
 408     temperatureCouplingImpl_->writeCheckpoint(
 409             checkpointData
 410                     ? std::make_optional(checkpointData->subCheckpointData("thermostat impl"))
 411                     : std::nullopt,
 412             cr);
 413 }
 414
 415 void VelocityScalingTemperatureCoupling::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
 416                                                                 const t_commrec* cr)
 417 {
 418     if (MASTER(cr))
 419     {
 420         doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
 421     }
 422     if (DOMAINDECOMP(cr))
 423     {
 424         dd_bcast(cr->dd, ssize(temperatureCouplingIntegral_) * int(sizeof(double)),
 425                  temperatureCouplingIntegral_.data());
 426     }
 427     temperatureCouplingImpl_->readCheckpoint(
 428             checkpointData
 429                     ? std::make_optional(checkpointData->subCheckpointData("thermostat impl"))
 430                     : std::nullopt,
 431             cr);
 432 }
 433
 434 const std::string& VelocityScalingTemperatureCoupling::clientID()
 435 {
 436     return identifier_;
 437 }
 438
 439 real VelocityScalingTemperatureCoupling::conservedEnergyContribution() const
 440 {
 441     if (reportPreviousConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
 442     {
 443         return std::accumulate(temperatureCouplingIntegralPreviousStep_.begin(),
 444                                temperatureCouplingIntegralPreviousStep_.end(), 0.0);
 445     }
 446     else
 447     {
 448         return std::accumulate(temperatureCouplingIntegral_.begin(),
 449                                temperatureCouplingIntegral_.end(), 0.0);
 450     }
 451 }
 452
 453 ISimulatorElement* VelocityScalingTemperatureCoupling::getElementPointerImpl(
 454         LegacySimulatorData*                    legacySimulatorData,
 455         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 456         StatePropagatorData gmx_unused* statePropagatorData,
 457         EnergyData gmx_unused*     energyData,
 458         FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
 459         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
 460         int                                   offset,
 461         UseFullStepKE                         useFullStepKE,
 462         ReportPreviousStepConservedEnergy     reportPreviousStepConservedEnergy)
 463 {
 464     // Element is now owned by the caller of this method, who will handle lifetime (see ModularSimulatorAlgorithm)
 465     auto* element = builderHelper->storeElement(std::make_unique<VelocityScalingTemperatureCoupling>(
 466             legacySimulatorData->inputrec->nsttcouple, offset, useFullStepKE, reportPreviousStepConservedEnergy,
 467             legacySimulatorData->inputrec->ld_seed, legacySimulatorData->inputrec->opts.ngtc,
 468             legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nsttcouple,
 469             legacySimulatorData->inputrec->opts.ref_t, legacySimulatorData->inputrec->opts.tau_t,
 470             legacySimulatorData->inputrec->opts.nrdf, energyData, legacySimulatorData->inputrec->etc));
 471     auto* thermostat = static_cast<VelocityScalingTemperatureCoupling*>(element);
 472     // Capturing pointer is safe because lifetime is handled by caller
 473     builderHelper->registerThermostat([thermostat](const PropagatorThermostatConnection& connection) {
 474         thermostat->connectWithPropagator(connection);
 475     });
 476     return element;
 477 }
 478
 479 } // namespace gmx