src/gromacs/modularsimulator/energydata.h

   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 Declares the energy element for the modular simulator
  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_ENERGYELEMENT_MICROSTATE_H
  45 #define GMX_ENERGYELEMENT_MICROSTATE_H
  46
  47 #include "gromacs/mdtypes/state.h"
  48
  49 #include "modularsimulatorinterfaces.h"
  50
  51 class gmx_ekindata_t;
  52 struct gmx_enerdata_t;
  53 struct gmx_mtop_t;
  54 struct ObservablesHistory;
  55 struct pull_t;
  56 struct t_fcdata;
  57 struct t_inputrec;
  58 struct SimulationGroups;
  59
  60 namespace gmx
  61 {
  62 enum class StartingBehavior;
  63 class Constraints;
  64 class EnergyOutput;
  65 class FreeEnergyPerturbationData;
  66 class GlobalCommunicationHelper;
  67 class LegacySimulatorData;
  68 class MDAtoms;
  69 class ModularSimulatorAlgorithmBuilderHelper;
  70 class ObservablesReducer;
  71 class ParrinelloRahmanBarostat;
  72 class StatePropagatorData;
  73 class VelocityScalingTemperatureCoupling;
  74 struct MDModulesNotifiers;
  75
  76 //! Function type for elements contributing energy
  77 using EnergyContribution = std::function<real(Step, Time)>;
  78
  79 /*! \internal
  80  * \ingroup module_modularsimulator
  81  * \brief Data class managing energies
  82  *
  83  * The EnergyData owns the EnergyObject,
  84  * the tensors for the different virials and the pressure as well as
  85  * the total dipole vector. It has a member class which is part of the
  86  * simulator loop and and is responsible
  87  * for saving energy data and writing it to trajectory.
  88  *
  89  * The EnergyData offers an interface to add virial contributions,
  90  * but also allows access to the raw pointers to tensor data, the
  91  * dipole vector, and the legacy energy data structures.
  92  *
  93  * The EnergyData owns an object of type EnergyData::Element,
  94  * which takes part in the simulation loop, allowing to record
  95  * and output energies during the simulation.
  96  */
  97 class EnergyData final
  98 {
  99 public:
 100     //! Constructor
 101     EnergyData(StatePropagatorData*        statePropagatorData,
 102                FreeEnergyPerturbationData* freeEnergyPerturbationData,
 103                const gmx_mtop_t&           globalTopology,
 104                const t_inputrec*           inputrec,
 105                const MDAtoms*              mdAtoms,
 106                gmx_enerdata_t*             enerd,
 107                gmx_ekindata_t*             ekind,
 108                const Constraints*          constr,
 109                FILE*                       fplog,
 110                t_fcdata*                   fcd,
 111                const MDModulesNotifiers&   mdModulesNotifiers,
 112                bool                        isMasterRank,
 113                ObservablesHistory*         observablesHistory,
 114                StartingBehavior            startingBehavior,
 115                bool                        simulationsShareState,
 116                pull_t*                     pullWork);
 117
 118     /*! \brief Final output
 119      *
 120      * Prints the averages to log. This is called from ModularSimulatorAlgorithm.
 121      *
 122      * \see ModularSimulatorAlgorithm::teardown
 123      */
 124     void teardown();
 125
 126     /*! \brief Add contribution to force virial
 127      *
 128      * This automatically resets the tensor if the step is higher
 129      * than the current step, starting the tensor calculation for
 130      * a new step at zero. Otherwise, it adds the new contribution
 131      * to the existing virial.
 132      */
 133     void addToForceVirial(const tensor virial, Step step);
 134
 135     /*! \brief Add contribution to constraint virial
 136      *
 137      * This automatically resets the tensor if the step is higher
 138      * than the current step, starting the tensor calculation for
 139      * a new step at zero. Otherwise, it adds the new contribution
 140      * to the existing virial.
 141      */
 142     void addToConstraintVirial(const tensor virial, Step step);
 143
 144     /*! \brief Get pointer to force virial tensor
 145      *
 146      * Allows access to the raw pointer to the tensor.
 147      */
 148     rvec* forceVirial(Step step);
 149
 150     /*! \brief Get pointer to constraint virial tensor
 151      *
 152      * Allows access to the raw pointer to the tensor.
 153      */
 154     rvec* constraintVirial(Step step);
 155
 156     /*! \brief Get pointer to total virial tensor
 157      *
 158      * Allows access to the raw pointer to the tensor.
 159      */
 160     rvec* totalVirial(Step step);
 161
 162     /*! \brief Get pointer to pressure tensor
 163      *
 164      * Allows access to the raw pointer to the tensor.
 165      */
 166     rvec* pressure(Step step);
 167
 168     /*! \brief Get pointer to mu_tot
 169      *
 170      * Allows access to the raw pointer to the dipole vector.
 171      */
 172     real* muTot();
 173
 174     /*! \brief Get pointer to energy structure
 175      *
 176      */
 177     gmx_enerdata_t* enerdata();
 178
 179     /*! \brief Get const pointer to energy structure
 180      *
 181      */
 182     const gmx_enerdata_t* enerdata() const;
 183
 184     /*! \brief Get pointer to kinetic energy structure
 185      *
 186      */
 187     gmx_ekindata_t* ekindata();
 188
 189     /*! \brief Get pointer to needToSumEkinhOld
 190      *
 191      */
 192     bool* needToSumEkinhOld();
 193
 194     /*! \brief Whether kinetic energy was read from checkpoint
 195      *
 196      * This is needed by the compute globals element
 197      * TODO: Remove this when moving global reduction to client system (#3421)
 198      */
 199     [[nodiscard]] bool hasReadEkinFromCheckpoint() const;
 200
 201     /*! \brief Add conserved energy contribution
 202      *
 203      * This allows other elements to register callbacks for contributions to
 204      * the conserved energy term.
 205      */
 206     void addConservedEnergyContribution(EnergyContribution&& energyContribution);
 207
 208     /*! \brief set Parrinello-Rahman barostat
 209      *
 210      * This allows to set a pointer to the Parrinello-Rahman barostat used to
 211      * print the box velocities.
 212      */
 213     void setParrinelloRahmanBoxVelocities(std::function<const rvec*()>&& parrinelloRahmanBoxVelocities);
 214
 215     /*! \brief Initialize energy history
 216      *
 217      * Kept as a static function to allow usage from legacy code
 218      * \todo Make member function once legacy use is not needed anymore
 219      */
 220     static void initializeEnergyHistory(StartingBehavior    startingBehavior,
 221                                         ObservablesHistory* observablesHistory,
 222                                         EnergyOutput*       energyOutput);
 223
 224     /*! \brief Request (local) kinetic energy update
 225      */
 226     void updateKineticEnergy();
 227
 228     //! The element taking part in the simulator loop
 229     class Element;
 230     //! Get pointer to element (whose lifetime is managed by this)
 231     Element* element();
 232
 233 private:
 234     /*! \brief Setup (needs file pointer)
 235      *
 236      * Initializes the EnergyOutput object, and does some logging output.
 237      *
 238      * \param mdoutf  File pointer
 239      */
 240     void setup(gmx_mdoutf* mdoutf);
 241
 242     /*! \brief Save data at energy steps
 243      *
 244      * \param step  The current step
 245      * \param time  The current time
 246      * \param isEnergyCalculationStep  Whether the current step is an energy calculation step
 247      * \param isFreeEnergyCalculationStep  Whether the current step is a free energy calculation step
 248      */
 249     void doStep(Step step, Time time, bool isEnergyCalculationStep, bool isFreeEnergyCalculationStep);
 250
 251     /*! \brief Write to energy trajectory
 252      *
 253      * This is only called by master - writes energy to trajectory and to log.
 254      */
 255     void write(gmx_mdoutf* outf, Step step, Time time, bool writeTrajectory, bool writeLog);
 256
 257     /*
 258      * Data owned by EnergyData
 259      */
 260     //! The element
 261     std::unique_ptr<Element> element_;
 262     //! The energy output object
 263     std::unique_ptr<EnergyOutput> energyOutput_;
 264     //! Helper object to checkpoint kinetic energy data
 265     ekinstate_t ekinstate_;
 266
 267     //! Whether this is the master rank
 268     const bool isMasterRank_;
 269
 270     //! The force virial tensor
 271     tensor forceVirial_;
 272     //! The constraint virial tensor
 273     tensor shakeVirial_;
 274     //! The total virial tensor
 275     tensor totalVirial_;
 276     //! The pressure tensor
 277     tensor pressure_;
 278     //! The total dipole moment
 279     rvec muTot_;
 280
 281     //! The step number of the current force virial tensor
 282     Step forceVirialStep_;
 283     //! The step number of the current constraint virial tensor
 284     Step shakeVirialStep_;
 285     //! The step number of the current total virial tensor
 286     Step totalVirialStep_;
 287     //! The step number of the current pressure tensor
 288     Step pressureStep_;
 289
 290     //! Whether ekinh_old needs to be summed up (set by compute globals)
 291     bool needToSumEkinhOld_;
 292     //! Whether we have read ekin from checkpoint
 293     bool hasReadEkinFromCheckpoint_;
 294
 295     //! Describes how the simulation (re)starts
 296     const StartingBehavior startingBehavior_;
 297
 298     /*
 299      * Pointers to Simulator data
 300      */
 301     // TODO: Clarify relationship to data objects and find a more robust alternative to raw pointers (#3583)
 302     //! Pointer to the state propagator data
 303     StatePropagatorData* statePropagatorData_;
 304     //! Pointer to the free energy perturbation data
 305     FreeEnergyPerturbationData* freeEnergyPerturbationData_;
 306
 307     //! Callbacks contributing to the conserved energy term
 308     std::vector<EnergyContribution> conservedEnergyContributions_;
 309     //! Callback to the Parrinello-Rahman box velocities
 310     std::function<const rvec*()> parrinelloRahmanBoxVelocities_;
 311
 312     //! Contains user input mdp options.
 313     const t_inputrec* inputrec_;
 314     //! Full system topology.
 315     const gmx_mtop_t& top_global_;
 316     //! Atom parameters for this domain.
 317     const MDAtoms* mdAtoms_;
 318     //! Energy data structure
 319     gmx_enerdata_t* enerd_;
 320     //! Kinetic energy data
 321     gmx_ekindata_t* ekind_;
 322     //! Handles constraints.
 323     const Constraints* constr_;
 324     //! Handles logging.
 325     FILE* fplog_;
 326     //! Helper struct for force calculations.
 327     t_fcdata* fcd_;
 328     //! Notifiers to MD modules
 329     const MDModulesNotifiers& mdModulesNotifiers_;
 330     //! Global topology groups
 331     const SimulationGroups* groups_;
 332     //! History of simulation observables.
 333     ObservablesHistory* observablesHistory_;
 334     //! Whether simulations share the state
 335     bool simulationsShareState_;
 336     //! The pull work object.
 337     pull_t* pullWork_;
 338 };
 339
 340 /*! \internal
 341  * \ingroup module_modularsimulator
 342  * \brief Element for EnergyData
 343  *
 344  * This member class allows EnergyData to take part in the simulator
 345  * loop.
 346  *
 347  * It subscribes to the trajectory signaller, the energy signaller,
 348  * and the logging signaller to know when an energy calculation is
 349  * needed and when a non-recording step is enough. The simulator
 350  * builder is responsible to place the element in a location at
 351  * which a valid energy state is available. The EnergyData::Element is
 352  * also a subscriber to the trajectory writer element, as it is
 353  * responsible to write energy data to trajectory.
 354  */
 355 class EnergyData::Element final :
 356     public ISimulatorElement,
 357     public ITrajectoryWriterClient,
 358     public ITrajectorySignallerClient,
 359     public IEnergySignallerClient,
 360     public ICheckpointHelperClient
 361 {
 362 public:
 363     //! Constructor
 364     Element(EnergyData* energyData, bool isMasterRank);
 365
 366     /*! \brief Register run function for step / time
 367      *
 368      * This needs to be called when the energies are at a full time step.
 369      * Positioning this element is the responsibility of the programmer.
 370      *
 371      * This is also the place at which the current state becomes the previous
 372      * state.
 373      *
 374      * \param step                 The step number
 375      * \param time                 The time
 376      * \param registerRunFunction  Function allowing to register a run function
 377      */
 378     void scheduleTask(Step step, Time time, const RegisterRunFunction& registerRunFunction) override;
 379
 380     //! No element setup needed
 381     void elementSetup() override {}
 382
 383     //! No element teardown needed
 384     void elementTeardown() override {}
 385
 386     //! ICheckpointHelperClient write checkpoint implementation
 387     void saveCheckpointState(std::optional<WriteCheckpointData> checkpointData, const t_commrec* cr) override;
 388     //! ICheckpointHelperClient read checkpoint implementation
 389     void restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData, const t_commrec* cr) override;
 390     //! ICheckpointHelperClient key implementation
 391     const std::string& clientID() override;
 392
 393     /*! \brief Factory method implementation
 394      *
 395      * \param legacySimulatorData  Pointer allowing access to simulator level data
 396      * \param builderHelper  ModularSimulatorAlgorithmBuilder helper object
 397      * \param statePropagatorData  Pointer to the \c StatePropagatorData object
 398      * \param energyData  Pointer to the \c EnergyData object
 399      * \param freeEnergyPerturbationData  Pointer to the \c FreeEnergyPerturbationData object
 400      * \param globalCommunicationHelper   Pointer to the \c GlobalCommunicationHelper object
 401      * \param observablesReducer          Pointer to the \c ObservablesReducer object
 402      *
 403      * \return  Pointer to the element to be added. Element needs to have been stored using \c storeElement
 404      */
 405     static ISimulatorElement* getElementPointerImpl(LegacySimulatorData* legacySimulatorData,
 406                                                     ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 407                                                     StatePropagatorData*        statePropagatorData,
 408                                                     EnergyData*                 energyData,
 409                                                     FreeEnergyPerturbationData* freeEnergyPerturbationData,
 410                                                     GlobalCommunicationHelper* globalCommunicationHelper,
 411                                                     ObservablesReducer*        observablesReducer);
 412
 413 private:
 414     EnergyData* energyData_;
 415
 416     /*! \brief Setup (needs file pointer)
 417      *
 418      * ITrajectoryWriterClient implementation.
 419      *
 420      * Initializes the EnergyOutput object, and does some logging output.
 421      *
 422      * \param mdoutf  File pointer
 423      */
 424     void trajectoryWriterSetup(gmx_mdoutf* mdoutf) override;
 425     //! No trajectory writer teardown needed
 426     void trajectoryWriterTeardown(gmx_mdoutf gmx_unused* outf) override {}
 427
 428     //! ITrajectoryWriterClient implementation.
 429     std::optional<SignallerCallback> registerTrajectorySignallerCallback(TrajectoryEvent event) override;
 430     //! ITrajectorySignallerClient implementation
 431     std::optional<ITrajectoryWriterCallback> registerTrajectoryWriterCallback(TrajectoryEvent event) override;
 432     //! IEnergySignallerClient implementation
 433     std::optional<SignallerCallback> registerEnergyCallback(EnergySignallerEvent event) override;
 434
 435
 436     //! CheckpointHelper identifier
 437     const std::string identifier_ = "EnergyElement";
 438     //! Helper function to read from / write to CheckpointData
 439     template<CheckpointDataOperation operation>
 440     void doCheckpointData(CheckpointData<operation>* checkpointData);
 441
 442     //! Whether this is the master rank
 443     const bool isMasterRank_;
 444     //! The next communicated energy writing step
 445     Step energyWritingStep_;
 446     //! The next communicated energy calculation step
 447     Step energyCalculationStep_;
 448     //! The next communicated free energy calculation step
 449     Step freeEnergyCalculationStep_;
 450 };
 451
 452 } // namespace gmx
 453
 454 #endif // GMX_ENERGYELEMENT_MICROSTATE_H