src/gromacs/modularsimulator/mttk.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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 Declares classes related to MTTK pressure coupling
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  *
  41  * This header is only used within the modular simulator module
  42  */
  43
  44 #ifndef GMX_MODULARSIMULATOR_MTTK_H
  45 #define GMX_MODULARSIMULATOR_MTTK_H
  46
  47 #include <queue>
  48
  49 #include "modularsimulatorinterfaces.h"
  50 #include "statepropagatordata.h"
  51
  52 namespace gmx
  53 {
  54 class EnergyData;
  55 class MttkPropagatorConnection;
  56 enum class ScheduleOnInitStep;
  57
  58 /*! \internal
  59  * \brief Struct collecting the propagator tags and offsets used by the MTTK pressure coupling
  60  *
  61  * MTTK scales the positions before and after the propagation step, and the
  62  * velocities before and after each propagation half step.
  63  */
  64 struct MttkPropagatorConnectionDetails
  65 {
  66     //! The tag of the position scaling before the propagation step
  67     PropagatorTag propagatorTagPrePosition;
  68     //! The tag of the position scaling after the propagation step
  69     PropagatorTag propagatorTagPostPosition;
  70     //! The offset at which position scaling is applied
  71     Offset positionOffset;
  72     //! The tag of the velocity scaling before the first propagation half step
  73     PropagatorTag propagatorTagPreVelocity1;
  74     //! The tag of the velocity scaling after the first propagation half step
  75     PropagatorTag propagatorTagPostVelocity1;
  76     //! The offset at which the first velocity half step scaling is applied
  77     Offset velocityOffset1;
  78     //! The tag of the velocity scaling before the second propagation half step
  79     PropagatorTag propagatorTagPreVelocity2;
  80     //! The tag of the velocity scaling after the second propagation half step
  81     PropagatorTag propagatorTagPostVelocity2;
  82     //! The offset at which the second velocity half step scaling is applied
  83     Offset velocityOffset2;
  84 };
  85
  86 /*! \internal
  87  * \brief Class holding the extra dof and parameters used by the MTTK algorithm
  88  *
  89  * As the Trotter update is split in several sub-steps (i.e. is updated
  90  * by several element instances), the MTTK degrees of freedom must be
  91  * stored centrally rather than by the single elements.
  92  *
  93  * This class manages these extra degrees of freedom. It controls access, keeps
  94  * track of the current time stamp of the dofs, calculates the energy
  95  * related to the dof at the requested times, and writes the data needed
  96  * for restarts to checkpoint. As this is not implementing the
  97  * ISimulatorElement interface, it is not part of the simulator loop, but
  98  * relies on callbacks to perform its duties.
  99  */
 100 class MttkData final : public ICheckpointHelperClient
 101 {
 102 public:
 103     //! Constructor
 104     MttkData(real                       referenceTemperature,
 105              real                       referencePressure,
 106              real                       couplingTimeStep,
 107              real                       couplingTime,
 108              real                       initialVolume,
 109              real                       numDegreesOfFreedom,
 110              real                       simulationTimeStep,
 111              const tensor               compressibility,
 112              const StatePropagatorData* statePropagatorData,
 113              MttkPropagatorConnection*  mttkPropagatorConnection);
 114
 115     //! Explicit copy constructor (interface has a standard destructor)
 116     MttkData(const MttkData& other);
 117
 118     //! The current kinetic energy of the MTTK degree of freedom
 119     [[nodiscard]] real kineticEnergy() const;
 120     /*! \brief Scale the MTTK dof velocity
 121      *
 122      * \param scalingFactor  The factor by which the velocity is scaled
 123      * \param scalingAtFullCouplingTimeStep  Whether the calling object is at full timestep
 124      *                                       (determines whether the integral is calculated)
 125      */
 126     void scale(real scalingFactor, bool scalingAtFullCouplingTimeStep);
 127
 128     //! The current MTTK dof velocity
 129     [[nodiscard]] real etaVelocity() const;
 130     //! Set a new MTTK velocity
 131     void setEtaVelocity(real etaVelocity, real etaVelocityTimeIncrement);
 132
 133     //! Get the inverse mass of the MTTK degree of freedom
 134     [[nodiscard]] real invEtaMass() const;
 135     //! Get the reference pressure
 136     real referencePressure() const;
 137     //! Pointer to the box velocities
 138     [[nodiscard]] rvec* boxVelocities();
 139
 140     //! Inform the propagators that scaling is needed
 141     void propagatorCallback(Step step) const;
 142
 143     //! ICheckpointHelperClient write checkpoint implementation
 144     void saveCheckpointState(std::optional<WriteCheckpointData> checkpointData, const t_commrec* cr) override;
 145     //! ICheckpointHelperClient read checkpoint implementation
 146     void restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData, const t_commrec* cr) override;
 147     //! ICheckpointHelperClient key implementation
 148     const std::string& clientID() override;
 149
 150     //! Build object and store in builder helper object
 151     static void build(LegacySimulatorData*                    legacySimulatorData,
 152                       ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 153                       StatePropagatorData*                    statePropagatorData,
 154                       EnergyData*                             energyData,
 155                       const MttkPropagatorConnectionDetails&  mttkPropagatorConnectionDetails);
 156
 157     //! Calculate the current value of the MTTK conserved energy if it is needed
 158     void calculateIntegralIfNeeded();
 159
 160     //! Identifier used to store objects
 161     static std::string dataID();
 162
 163 private:
 164     //! Return the current value of the MTTK dof contribution to the conserved energy
 165     double temperatureCouplingIntegral(Time time) const;
 166     //! Update the position and velocity scaling factors
 167     void updateScalingFactors();
 168     //! Update the reference temperature
 169     void updateReferenceTemperature(real temperature, ReferenceTemperatureChangeAlgorithm algorithm);
 170     //! Calculate the MTTK conserved energy
 171     void calculateIntegral(real volume);
 172
 173     //! The coupling time step
 174     const real couplingTimeStep_;
 175     //! The velocity of the MTTK dof
 176     real etaVelocity_;
 177     //! The inverse mass of the MTTK dof
 178     real invMass_;
 179     //! The current time stamp of the MTTK dof velocity
 180     Time etaVelocityTime_;
 181     //! The current value of the MTTK dof contribution to the conserved energy
 182     double temperatureCouplingIntegral_;
 183     //! The current time stamp of the conserved energy contribution
 184     Time integralTime_;
 185     //! The reference pressure
 186     const real referencePressure_;
 187     //! The current box velocities (used for reporting only)
 188     tensor boxVelocity_;
 189     //! The number of degrees of freedom in the first temperature group
 190     const real numDegreesOfFreedom_;
 191     //! The simulation time step - by how much the propagators are advancing the positions
 192     const real simulationTimeStep_;
 193     //! The reference temperature the mass is based on
 194     real referenceTemperature_;
 195
 196     // TODO: Clarify relationship to data objects and find a more robust alternative to raw pointers (#3583)
 197     //! Pointer to the micro state data (access to the current box)
 198     const StatePropagatorData* statePropagatorData_;
 199     //! Pointer to the propagator connection object
 200     MttkPropagatorConnection* mttkPropagatorConnection_;
 201
 202     //! CheckpointHelper identifier
 203     const std::string identifier_ = "MttkData";
 204     //! Helper function to read from / write to CheckpointData
 205     template<CheckpointDataOperation operation>
 206     void doCheckpointData(CheckpointData<operation>* checkpointData);
 207 };
 208
 209 /*! \internal
 210  * \brief Object holding the callbacks and scaling views for the connection
 211  *        of MTTKElement objects to propagators
 212  *
 213  * As the Trotter update is split in several sub-steps (i.e. is updated
 214  * by several element instances), the connection to propagators must be
 215  * stored centrally rather than by the single elements.
 216  */
 217 class MttkPropagatorConnection
 218 {
 219 public:
 220     //! Build object and store in builder helper object
 221     static void build(ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 222                       const PropagatorTag&                    propagatorTagPrePosition,
 223                       const PropagatorTag&                    propagatorTagPostPosition,
 224                       int                                     positionOffset,
 225                       const PropagatorTag&                    propagatorTagPreVelocity1,
 226                       const PropagatorTag&                    propagatorTagPostVelocity1,
 227                       int                                     velocityOffset1,
 228                       const PropagatorTag&                    propagatorTagPreVelocity2,
 229                       const PropagatorTag&                    propagatorTagPostVelocity2,
 230                       int                                     velocityOffset2);
 231
 232     //! Call the propagator callbacks
 233     void propagatorCallback(Step step) const;
 234
 235     //! Set position scaling
 236     void setPositionScaling(real preStepScaling, real postStepScaling);
 237     //! Set velocity scaling
 238     void setVelocityScaling(real preStepScaling, real postStepScaling);
 239
 240     //! Identifier used to store objects
 241     static std::string dataID();
 242
 243 private:
 244     //! Connect to pre-step velocity scaling
 245     void connectWithPropagatorVelocityPreStepScaling(const PropagatorConnection& connectionData,
 246                                                      const PropagatorTag&        propagatorTag,
 247                                                      int                         offset);
 248     //! Connect to post-step velocity scaling
 249     void connectWithPropagatorVelocityPostStepScaling(const PropagatorConnection& connectionData,
 250                                                       const PropagatorTag&        propagatorTag,
 251                                                       int                         offset);
 252     //! Connect to pre-step position scaling
 253     void connectWithPropagatorPositionPreStepScaling(const PropagatorConnection& connectionData,
 254                                                      const PropagatorTag&        propagatorTag,
 255                                                      int                         offset);
 256     //! Connect to post-step position scaling
 257     void connectWithPropagatorPositionPostStepScaling(const PropagatorConnection& connectionData,
 258                                                       const PropagatorTag&        propagatorTag,
 259                                                       int                         offset);
 260
 261     //! View on the scaling factors of the position propagators (before step)
 262     std::vector<ArrayRef<real>> startPositionScalingFactors_;
 263     //! View on the scaling factors of the position propagators (after step)
 264     std::vector<ArrayRef<real>> endPositionScalingFactors_;
 265     //! View on the scaling factors of the velocity propagators (before step)
 266     std::vector<ArrayRef<real>> startVelocityScalingFactors_;
 267     //! View on the scaling factors of the velocity propagators (after step)
 268     std::vector<ArrayRef<real>> endVelocityScalingFactors_;
 269     //! Callbacks to let propagators know that we will update scaling factor, plus their offset
 270     std::vector<std::tuple<PropagatorCallback, int>> propagatorCallbacks_;
 271 };
 272
 273 /*! \internal
 274  * \brief Element propagating the MTTK degree of freedom
 275  *
 276  * This roughly follows Tuckerman et al. 2006 (DOI: 10.1088/0305-4470/39/19/S18)
 277  *
 278  * This currently only implements the isotropic version of MTTK and does
 279  * not work with constraints. It also adopts some conventions chosen by
 280  * the legacy GROMACS implementation, such as using the number of degrees
 281  * of freedom and the kinetic energy scaling of the first temperature group
 282  * to calculate the current pressure. (It does use the full system kinetic
 283  * energy, however!)
 284  */
 285 class MttkElement final : public ISimulatorElement
 286 {
 287 public:
 288     //! Constructor
 289     MttkElement(int                        nstcouple,
 290                 int                        offset,
 291                 real                       propagationTimeStep,
 292                 ScheduleOnInitStep         scheduleOnInitStep,
 293                 Step                       initStep,
 294                 const StatePropagatorData* statePropagatorData,
 295                 EnergyData*                energyData,
 296                 MttkData*                  mttkData,
 297                 PbcType                    pbcType,
 298                 int                        numWalls,
 299                 real                       numDegreesOfFreedom);
 300
 301     /*! \brief Register run function for step / time
 302      *
 303      * \param step                 The step number
 304      * \param time                 The time
 305      * \param registerRunFunction  Function allowing to register a run function
 306      */
 307     void scheduleTask(Step step, Time time, const RegisterRunFunction& registerRunFunction) override;
 308
 309     //! No element setup needed
 310     void elementSetup() override {}
 311     //! No element teardown needed
 312     void elementTeardown() override {}
 313
 314     /*! \brief Factory method implementation
 315      *
 316      * \param legacySimulatorData  Pointer allowing access to simulator level data
 317      * \param builderHelper  ModularSimulatorAlgorithmBuilder helper object
 318      * \param statePropagatorData  Pointer to the \c StatePropagatorData object
 319      * \param energyData  Pointer to the \c EnergyData object
 320      * \param freeEnergyPerturbationData  Pointer to the \c FreeEnergyPerturbationData object
 321      * \param globalCommunicationHelper  Pointer to the \c GlobalCommunicationHelper object
 322      * \param observablesReducer  Pointer to the \c ObservablesReducer object
 323      * \param offset  The step offset at which the thermostat is applied
 324      * \param scheduleOnInitStep  Whether the element is scheduled on the initial step
 325      * \param mttkPropagatorConnectionDetails  Reference to the \c MttkPropagatorConnectionDetails object containing propagator tags and offsets
 326      *
 327      * \return  Pointer to the element to be added. Element needs to have been stored using \c storeElement
 328      */
 329     static ISimulatorElement*
 330     getElementPointerImpl(LegacySimulatorData*                    legacySimulatorData,
 331                           ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 332                           StatePropagatorData*                    statePropagatorData,
 333                           EnergyData*                             energyData,
 334                           FreeEnergyPerturbationData*             freeEnergyPerturbationData,
 335                           GlobalCommunicationHelper*              globalCommunicationHelper,
 336                           ObservablesReducer*                     observablesReducer,
 337                           Offset                                  offset,
 338                           ScheduleOnInitStep                      scheduleOnInitStep,
 339                           const MttkPropagatorConnectionDetails&  mttkPropagatorConnectionDetails);
 340
 341 private:
 342     //! Propagation of the MTTK dof velocity
 343     void propagateEtaVelocity(Step step);
 344
 345     //! The periodic boundary condition type
 346     const PbcType pbcType_;
 347     //! The number of walls
 348     const int numWalls_;
 349     //! The number of degrees of freedom in the first temperature group
 350     const real numDegreesOfFreedom_;
 351
 352     //! The period at which the thermostat is applied
 353     const int nstcouple_;
 354     //! If != 0, offset the step at which the thermostat is applied
 355     const int offset_;
 356     //! The propagation time step - by how much we propagate the MTTK dof
 357     const real propagationTimeStep_;
 358     //! Whether we're scheduling on the first step
 359     const ScheduleOnInitStep scheduleOnInitStep_;
 360     //! The initial step number
 361     const Step initialStep_;
 362
 363     // TODO: Clarify relationship to data objects and find a more robust alternative to raw pointers (#3583)
 364     //! Pointer to the micro state
 365     const StatePropagatorData* statePropagatorData_;
 366     //! Pointer to the energy data (for ekindata)
 367     EnergyData* energyData_;
 368     //! Pointer to the MTTK data
 369     MttkData* mttkData_;
 370 };
 371
 372 /*! \internal
 373  * \brief This element scales the box based on the MTTK dof
 374  */
 375 class MttkBoxScaling final : public ISimulatorElement
 376 {
 377 public:
 378     //! Constructor
 379     MttkBoxScaling(real simulationTimeStep, StatePropagatorData* statePropagatorData, MttkData* mttkData);
 380
 381     /*! \brief Register run function for step / time
 382      *
 383      * \param step                 The step number
 384      * \param time                 The time
 385      * \param registerRunFunction  Function allowing to register a run function
 386      */
 387     void scheduleTask(Step step, Time time, const RegisterRunFunction& registerRunFunction) override;
 388
 389     //! Sanity check at setup time
 390     void elementSetup() override {}
 391     //! No element teardown needed
 392     void elementTeardown() override {}
 393
 394     /*! \brief Factory method implementation
 395      *
 396      * \param legacySimulatorData  Pointer allowing access to simulator level data
 397      * \param builderHelper  ModularSimulatorAlgorithmBuilder helper object
 398      * \param statePropagatorData  Pointer to the \c StatePropagatorData object
 399      * \param energyData  Pointer to the \c EnergyData object
 400      * \param freeEnergyPerturbationData  Pointer to the \c FreeEnergyPerturbationData object
 401      * \param globalCommunicationHelper   Pointer to the \c GlobalCommunicationHelper object
 402      * \param observablesReducer          Pointer to the \c ObservablesReducer object
 403      * \param mttkPropagatorConnectionDetails  Reference to the \c MttkPropagatorConnectionDetails object containing propagator tags and offsets
 404      *
 405      * \return  Pointer to the element to be added. Element needs to have been stored using \c storeElement
 406      */
 407     static ISimulatorElement*
 408     getElementPointerImpl(LegacySimulatorData*                    legacySimulatorData,
 409                           ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 410                           StatePropagatorData*                    statePropagatorData,
 411                           EnergyData*                             energyData,
 412                           FreeEnergyPerturbationData*             freeEnergyPerturbationData,
 413                           GlobalCommunicationHelper*              globalCommunicationHelper,
 414                           ObservablesReducer*                     observablesReducer,
 415                           const MttkPropagatorConnectionDetails&  mttkPropagatorConnectionDetails);
 416
 417 private:
 418     //! Scale the box
 419     void scaleBox();
 420     //! The simulation time step
 421     const real simulationTimeStep_;
 422
 423     // TODO: Clarify relationship to data objects and find a more robust alternative to raw pointers (#3583)
 424     //! Pointer to the micro state
 425     StatePropagatorData* statePropagatorData_;
 426     //! Pointer to the MTTK data (nullptr if this is not connected to barostat)
 427     MttkData* mttkData_;
 428 };
 429
 430 } // namespace gmx
 431
 432 #endif // GMX_MODULARSIMULATOR_MTTK_H