src/gromacs/modularsimulator/energyelement.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, 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 /*! \libinternal \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
  42 #ifndef GMX_ENERGYELEMENT_MICROSTATE_H
  43 #define GMX_ENERGYELEMENT_MICROSTATE_H
  44
  45 #include "gromacs/math/vectypes.h"
  46 #include "gromacs/mdtypes/state.h"
  47
  48 #include "modularsimulatorinterfaces.h"
  49
  50 struct gmx_ekindata_t;
  51 struct gmx_enerdata_t;
  52 struct gmx_mtop_t;
  53 struct ObservablesHistory;
  54 struct t_fcdata;
  55 struct t_inputrec;
  56 struct SimulationGroups;
  57
  58 namespace gmx
  59 {
  60 enum class StartingBehavior;
  61 class Constraints;
  62 class EnergyOutput;
  63 class FreeEnergyPerturbationElement;
  64 class MDAtoms;
  65 class ParrinelloRahmanBarostat;
  66 class StatePropagatorData;
  67 class VRescaleThermostat;
  68 struct MdModulesNotifier;
  69
  70 /*! \libinternal
  71  * \ingroup module_modularsimulator
  72  * \brief Element managing energies
  73  *
  74  * The EnergyElement owns the EnergyObject, and is hence responsible
  75  * for saving energy data and writing it to trajectory. It also owns
  76  * the tensors for the different virials and the pressure as well as
  77  * the total dipole vector.
  78  *
  79  * It subscribes to the trajectory signaller, the energy signaller,
  80  * and the logging signaller to know when an energy calculation is
  81  * needed and when a non-recording step is enough. The simulator
  82  * builder is responsible to place the element in a location at
  83  * which a valid energy state is available. The EnergyElement is
  84  * also a subscriber to the trajectory writer element, as it is
  85  * responsible to write energy data to trajectory.
  86  *
  87  * The EnergyElement offers an interface to add virial contributions,
  88  * but also allows access to the raw pointers to tensor data, the
  89  * dipole vector, and the legacy energy data structures.
  90  */
  91 class EnergyElement final :
  92     public ISimulatorElement,
  93     public ITrajectoryWriterClient,
  94     public ITrajectorySignallerClient,
  95     public IEnergySignallerClient,
  96     public ICheckpointHelperClient
  97 {
  98 public:
  99     //! Constructor
 100     EnergyElement(StatePropagatorData*           statePropagatorData,
 101                   FreeEnergyPerturbationElement* freeEnergyPerturbationElement,
 102                   const gmx_mtop_t*              globalTopology,
 103                   const t_inputrec*              inputrec,
 104                   const MDAtoms*                 mdAtoms,
 105                   gmx_enerdata_t*                enerd,
 106                   gmx_ekindata_t*                ekind,
 107                   const Constraints*             constr,
 108                   FILE*                          fplog,
 109                   t_fcdata*                      fcd,
 110                   const MdModulesNotifier&       mdModulesNotifier,
 111                   bool                           isMasterRank,
 112                   ObservablesHistory*            observablesHistory,
 113                   StartingBehavior               startingBehavior);
 114
 115     /*! \brief Register run function for step / time
 116      *
 117      * This needs to be called when the energies are at a full time step.
 118      * Positioning this element is the responsibility of the programmer.
 119      *
 120      * This is also the place at which the current state becomes the previous
 121      * state.
 122      *
 123      * @param step                 The step number
 124      * @param time                 The time
 125      * @param registerRunFunction  Function allowing to register a run function
 126      */
 127     void scheduleTask(Step step, Time time, const RegisterRunFunctionPtr& registerRunFunction) override;
 128
 129     //! No element setup needed
 130     void elementSetup() override {}
 131
 132     /*! \brief Final output
 133      *
 134      * Prints the averages to log.
 135      */
 136     void elementTeardown() override;
 137
 138     /*! \brief Add contribution to force virial
 139      *
 140      * This automatically resets the tensor if the step is higher
 141      * than the current step, starting the tensor calculation for
 142      * a new step at zero. Otherwise, it adds the new contribution
 143      * to the existing virial.
 144      */
 145     void addToForceVirial(const tensor virial, Step step);
 146
 147     /*! \brief Add contribution to constraint virial
 148      *
 149      * This automatically resets the tensor if the step is higher
 150      * than the current step, starting the tensor calculation for
 151      * a new step at zero. Otherwise, it adds the new contribution
 152      * to the existing virial.
 153      */
 154     void addToConstraintVirial(const tensor virial, Step step);
 155
 156     /*! \brief Get pointer to force virial tensor
 157      *
 158      * Allows access to the raw pointer to the tensor.
 159      */
 160     rvec* forceVirial(Step step);
 161
 162     /*! \brief Get pointer to constraint virial tensor
 163      *
 164      * Allows access to the raw pointer to the tensor.
 165      */
 166     rvec* constraintVirial(Step step);
 167
 168     /*! \brief Get pointer to total virial tensor
 169      *
 170      * Allows access to the raw pointer to the tensor.
 171      */
 172     rvec* totalVirial(Step step);
 173
 174     /*! \brief Get pointer to pressure tensor
 175      *
 176      * Allows access to the raw pointer to the tensor.
 177      */
 178     rvec* pressure(Step step);
 179
 180     /*! \brief Get pointer to mu_tot
 181      *
 182      * Allows access to the raw pointer to the dipole vector.
 183      */
 184     real* muTot();
 185
 186     /*! \brief Get pointer to energy structure
 187      *
 188      */
 189     gmx_enerdata_t* enerdata();
 190
 191     /*! \brief Get pointer to kinetic energy structure
 192      *
 193      */
 194     gmx_ekindata_t* ekindata();
 195
 196     /*! \brief Get pointer to needToSumEkinhOld
 197      *
 198      */
 199     bool* needToSumEkinhOld();
 200
 201     /*! \brief set vrescale thermostat
 202      *
 203      * This allows to set a pointer to the vrescale thermostat used to
 204      * print the thermostat integral.
 205      * TODO: This should be made obsolete my a more modular energy element
 206      */
 207     void setVRescaleThermostat(const VRescaleThermostat* vRescaleThermostat);
 208
 209     /*! \brief set Parrinello-Rahman barostat
 210      *
 211      * This allows to set a pointer to the Parrinello-Rahman barostat used to
 212      * print the box velocities.
 213      * TODO: This should be made obsolete my a more modular energy element
 214      */
 215     void setParrinelloRahamnBarostat(const ParrinelloRahmanBarostat* parrinelloRahmanBarostat);
 216
 217     /*! \brief Initialize energy history
 218      *
 219      * Kept as a static function to allow usage from legacy code
 220      * \todo Make member function once legacy use is not needed anymore
 221      */
 222     static void initializeEnergyHistory(StartingBehavior    startingBehavior,
 223                                         ObservablesHistory* observablesHistory,
 224                                         EnergyOutput*       energyOutput);
 225
 226 private:
 227     /*! \brief Setup (needs file pointer)
 228      *
 229      * ITrajectoryWriterClient implementation.
 230      *
 231      * Initializes the EnergyOutput object, and does some logging output.
 232      *
 233      * @param mdoutf  File pointer
 234      */
 235     void trajectoryWriterSetup(gmx_mdoutf* mdoutf) override;
 236     //! No trajectory writer teardown needed
 237     void trajectoryWriterTeardown(gmx_mdoutf gmx_unused* outf) override {}
 238
 239     //! ITrajectoryWriterClient implementation.
 240     SignallerCallbackPtr registerTrajectorySignallerCallback(TrajectoryEvent event) override;
 241     //! ITrajectorySignallerClient implementation
 242     ITrajectoryWriterCallbackPtr registerTrajectoryWriterCallback(TrajectoryEvent event) override;
 243     //! IEnergySignallerClient implementation
 244     SignallerCallbackPtr registerEnergyCallback(EnergySignallerEvent event) override;
 245
 246     /*! \brief Save data at energy steps
 247      *
 248      * @param time  The current time
 249      * @param isEnergyCalculationStep  Whether the current step is an energy calculation step
 250      * @param isFreeEnergyCalculationStep  Whether the current step is a free energy calculation step
 251      */
 252     void doStep(Time time, bool isEnergyCalculationStep, bool isFreeEnergyCalculationStep);
 253
 254     /*! \brief Write to energy trajectory
 255      *
 256      * This is only called by master - writes energy to trajectory and to log.
 257      */
 258     void write(gmx_mdoutf* outf, Step step, Time time, bool writeTrajectory, bool writeLog);
 259
 260     //! ICheckpointHelperClient implementation
 261     void writeCheckpoint(t_state* localState, t_state* globalState) override;
 262
 263     /*
 264      * Data owned by EnergyElement
 265      */
 266     //! The energy output object
 267     std::unique_ptr<EnergyOutput> energyOutput_;
 268
 269     //! Whether this is the master rank
 270     const bool isMasterRank_;
 271     //! The next communicated energy writing step
 272     Step energyWritingStep_;
 273     //! The next communicated energy calculation step
 274     Step energyCalculationStep_;
 275     //! The next communicated free energy calculation step
 276     Step freeEnergyCalculationStep_;
 277
 278     //! The force virial tensor
 279     tensor forceVirial_;
 280     //! The constraint virial tensor
 281     tensor shakeVirial_;
 282     //! The total virial tensor
 283     tensor totalVirial_;
 284     //! The pressure tensor
 285     tensor pressure_;
 286     //! The total dipole moment
 287     rvec muTot_;
 288
 289     //! The step number of the current force virial tensor
 290     Step forceVirialStep_;
 291     //! The step number of the current constraint virial tensor
 292     Step shakeVirialStep_;
 293     //! The step number of the current total virial tensor
 294     Step totalVirialStep_;
 295     //! The step number of the current pressure tensor
 296     Step pressureStep_;
 297
 298     //! Whether ekinh_old needs to be summed up (set by compute globals)
 299     bool needToSumEkinhOld_;
 300
 301     //! Describes how the simulation (re)starts
 302     const StartingBehavior startingBehavior_;
 303
 304     //! Legacy state object used to communicate with energy output
 305     t_state dummyLegacyState_;
 306
 307     /*
 308      * Pointers to Simulator data
 309      */
 310     //! Pointer to the state propagator data
 311     StatePropagatorData* statePropagatorData_;
 312     //! Pointer to the free energy perturbation element
 313     FreeEnergyPerturbationElement* freeEnergyPerturbationElement_;
 314     //! Pointer to the vrescale thermostat
 315     const VRescaleThermostat* vRescaleThermostat_;
 316     //! Pointer to the Parrinello-Rahman barostat
 317     const ParrinelloRahmanBarostat* parrinelloRahmanBarostat_;
 318     //! Contains user input mdp options.
 319     const t_inputrec* inputrec_;
 320     //! Full system topology.
 321     const gmx_mtop_t* top_global_;
 322     //! Atom parameters for this domain.
 323     const MDAtoms* mdAtoms_;
 324     //! Energy data structure
 325     gmx_enerdata_t* enerd_;
 326     //! Kinetic energy data
 327     gmx_ekindata_t* ekind_;
 328     //! Handles constraints.
 329     const Constraints* constr_;
 330     //! Handles logging.
 331     FILE* fplog_;
 332     //! Helper struct for force calculations.
 333     t_fcdata* fcd_;
 334     //! Notification to MD modules
 335     const MdModulesNotifier& mdModulesNotifier_;
 336     //! Global topology groups
 337     const SimulationGroups* groups_;
 338     //! History of simulation observables.
 339     ObservablesHistory* observablesHistory_;
 340 };
 341
 342 } // namespace gmx
 343
 344 #endif // GMX_ENERGYELEMENT_MICROSTATE_H