src/gromacs/mdlib/energyoutput.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 /*! \libinternal \file
  36  * \brief Header for the code that writes energy-like quantities.
  37  *
  38  * \author Mark Abraham <mark.j.abraham@gmail.com>
  39  * \author Paul Bauer <paul.bauer.q@gmail.com>
  40  * \author Artem Zhmurov <zhmurov@gmail.com>
  41  *
  42  * \ingroup module_mdlib
  43  * \inlibraryapi
  44  */
  45 #ifndef GMX_MDLIB_ENERGYOUTPUT_H
  46 #define GMX_MDLIB_ENERGYOUTPUT_H
  47
  48 #include <cstdio>
  49
  50 #include <memory>
  51
  52 #include "gromacs/math/vectypes.h"
  53 #include "gromacs/mdtypes/md_enums.h"
  54 #include "gromacs/mdtypes/enerdata.h"
  55 #include "gromacs/topology/ifunc.h"
  56 #include "gromacs/utility/arrayref.h"
  57 #include "gromacs/utility/real.h"
  58
  59 class energyhistory_t;
  60 struct ener_file;
  61 class gmx_ekindata_t;
  62 struct gmx_enerdata_t;
  63 struct SimulationGroups;
  64 struct gmx_mtop_t;
  65 struct gmx_output_env_t;
  66 struct pull_t;
  67 struct t_ebin;
  68 struct t_expanded;
  69 struct t_fcdata;
  70 struct t_grpopts;
  71 struct t_inputrec;
  72 struct t_lambda;
  73 class t_state;
  74
  75 struct t_mde_delta_h_coll;
  76
  77 namespace gmx
  78 {
  79 class Awh;
  80 class Constraints;
  81 struct MDModulesNotifiers;
  82 enum class StartingBehavior;
  83 } // namespace gmx
  84
  85 //! \brief Printed names for intergroup energies
  86 const char* enumValueToString(NonBondedEnergyTerms enumValue);
  87
  88 /* \brief delta_h block type enum: the kinds of energies written out. */
  89 enum
  90 {
  91     //! Delta H BAR energy difference
  92     dhbtDH = 0,
  93     //! dH/dlambda derivative
  94     dhbtDHDL = 1,
  95     //! System energy
  96     dhbtEN,
  97     //! pV term
  98     dhbtPV,
  99     //! Expanded ensemble statistics
 100     dhbtEXPANDED,
 101     //! Total number of energy types in this enum
 102     dhbtNR
 103 };
 104
 105 namespace gmx
 106 {
 107 class EnergyDriftTracker;
 108
 109 /*! \internal
 110  * \brief Arrays connected to Pressure and Temperature coupling
 111  */
 112 struct PTCouplingArrays
 113 {
 114     //! Box velocities for Parrinello-Rahman P-coupling.
 115     const rvec* boxv;
 116     //! Nose-Hoover coordinates.
 117     ArrayRef<const double> nosehoover_xi;
 118     //! Nose-Hoover velocities.
 119     ArrayRef<const double> nosehoover_vxi;
 120     //! Pressure Nose-Hoover coordinates.
 121     ArrayRef<const double> nhpres_xi;
 122     //! Pressure Nose-Hoover velocities.
 123     ArrayRef<const double> nhpres_vxi;
 124 };
 125
 126 /* Energy output class
 127  *
 128  * This is the collection of energy averages collected during mdrun, and to
 129  * be written out to the .edr file.
 130  *
 131  * \todo Use more std containers.
 132  * \todo Remove GMX_CONSTRAINTVIR
 133  * \todo Write free-energy output also to energy file (after adding more tests)
 134  */
 135 class EnergyOutput
 136 {
 137 public:
 138     /*! \brief Initiate MD energy bin
 139      *
 140      * \param[in] fp_ene     Energy output file.
 141      * \param[in] mtop       Topology.
 142      * \param[in] inputrec   Input parameters.
 143      * \param[in] pull_work  Pulling simulations data
 144      * \param[in] fp_dhdl    FEP file.
 145      * \param[in] isRerun    Is this is a rerun instead of the simulations.
 146      * \param[in] startingBehavior  Run starting behavior.
 147      * \param[in] simulationsShareState  Tells whether the physical state is shared over simulations
 148      * \param[in] mdModulesNotifiers Notifications to MD modules.
 149      */
 150     EnergyOutput(ener_file*                fp_ene,
 151                  const gmx_mtop_t&         mtop,
 152                  const t_inputrec&         inputrec,
 153                  const pull_t*             pull_work,
 154                  FILE*                     fp_dhdl,
 155                  bool                      isRerun,
 156                  StartingBehavior          startingBehavior,
 157                  bool                      simulationsShareState,
 158                  const MDModulesNotifiers& mdModulesNotifiers);
 159
 160     ~EnergyOutput();
 161
 162     /*! \brief Update the averaging structures.
 163      *
 164      * Called every step on which the thermodynamic values are evaluated.
 165      *
 166      * \param[in] bDoDHDL           Whether the FEP is enabled.
 167      * \param[in] bSum              If this stepshould be recorded to compute sums and averages.
 168      * \param[in] time              Current simulation time.
 169      * \param[in] tmass             Total mass
 170      * \param[in] enerd             Energy data object.
 171      * \param[in] fep               FEP data.
 172      * \param[in] expand            Expanded ensemble (for FEP).
 173      * \param[in] lastbox           PBC data.
 174      * \param[in] ptCouplingArrays  Arrays connected to pressure and temperature coupling.
 175      * \param[in] fep_state         The current alchemical state we are in.
 176      * \param[in] svir              Constraint virial.
 177      * \param[in] fvir              Force virial.
 178      * \param[in] vir               Total virial.
 179      * \param[in] pres              Pressure.
 180      * \param[in] ekind             Kinetic energy data.
 181      * \param[in] mu_tot            Total dipole.
 182      * \param[in] constr            Constraints object to get RMSD from (for LINCS only).
 183      */
 184     void addDataAtEnergyStep(bool                    bDoDHDL,
 185                              bool                    bSum,
 186                              double                  time,
 187                              real                    tmass,
 188                              const gmx_enerdata_t*   enerd,
 189                              const t_lambda*         fep,
 190                              const t_expanded*       expand,
 191                              const matrix            lastbox,
 192                              PTCouplingArrays        ptCouplingArrays,
 193                              int                     fep_state,
 194                              const tensor            svir,
 195                              const tensor            fvir,
 196                              const tensor            vir,
 197                              const tensor            pres,
 198                              const gmx_ekindata_t*   ekind,
 199                              const rvec              mu_tot,
 200                              const gmx::Constraints* constr);
 201
 202     /*! \brief Update the data averaging structure counts.
 203      *
 204      * Updates the number of steps, the values have not being computed.
 205      */
 206     void recordNonEnergyStep();
 207
 208     /*! \brief Writes current quantites to log and energy files.
 209      *
 210      * Prints current values of energies, pressure, temperature, restraint
 211      * data, etc. to energy output file and to the log file (if not nullptr).
 212      *
 213      * This function only does something useful when bEne || bDR || bOR || log.
 214      *
 215      * \todo Perhaps this responsibility should involve some other
 216      *       object visiting all the contributing objects.
 217      *
 218      * \param[in] fp_ene   Energy file for the output.
 219      * \param[in] bEne     If it is a step for energy output or last step.
 220      * \param[in] bDR      If it is a step of writing distance restraints.
 221      * \param[in] bOR      If it is a step of writing orientation restraints.
 222      * \param[in] log      Pointer to the log file.
 223      * \param[in] step     Current step.
 224      * \param[in] time     Current simulation time.
 225      * \param[in] fcd      Bonded force computation data,
 226      *                     including orientation and distance restraints.
 227      * \param[in] awh      AWH data.
 228      */
 229     void printStepToEnergyFile(ener_file* fp_ene,
 230                                bool       bEne,
 231                                bool       bDR,
 232                                bool       bOR,
 233                                FILE*      log,
 234                                int64_t    step,
 235                                double     time,
 236                                t_fcdata*  fcd,
 237                                gmx::Awh*  awh);
 238
 239     /*! \brief Print reference temperatures for annealing groups.
 240      *
 241      * Nothing is done if log is nullptr.
 242      *
 243      * \param[in] log     Log file to print to.
 244      * \param[in] groups  Information on atom groups.
 245      * \param[in] opts    Atom temperature coupling groups options
 246      *                    (annealing is done by groups).
 247      */
 248     static void printAnnealingTemperatures(FILE* log, const SimulationGroups* groups, const t_grpopts* opts);
 249
 250     /*! \brief Prints average values to log file.
 251      *
 252      * This is called at the end of the simulation run to print accumulated average values.
 253      * Nothing it done if log is nullptr.
 254      *
 255      * \param[in]   log      Where to print.
 256      * \param[in]   groups   Atom groups.
 257      */
 258     void printAverages(FILE* log, const SimulationGroups* groups);
 259
 260     /*! \brief Get the number of thermodynamic terms recorded.
 261      *
 262      * \todo Refactor this to return the expected output size,
 263      *       rather than exposing the implementation details about
 264      *       thermodynamic terms.
 265      *
 266      * \returns Number of thermodynamic terms.
 267      */
 268     int numEnergyTerms() const;
 269
 270     /*! \brief Fill the energyhistory_t data.
 271      *
 272      * Between .edr writes, the averages are history dependent,
 273      * and that history needs to be retained in checkpoints.
 274      * These functions set/read the energyhistory_t class
 275      * that is written to checkpoints.
 276      *
 277      * \param[out] enerhist  Energy history data structure.
 278      */
 279     void fillEnergyHistory(energyhistory_t* enerhist) const;
 280
 281     /*! \brief Restore from energyhistory_t data.
 282      *
 283      * \param[in] enerhist  Energy history data structure.
 284      */
 285     void restoreFromEnergyHistory(const energyhistory_t& enerhist);
 286
 287     //! Print an output header to the log file.
 288     static void printHeader(FILE* log, int64_t steps, double time);
 289
 290     /*! \brief Print conserved energy drift message to \p fplog
 291      *
 292      * Note that this is only over the current run (times are printed),
 293      * this is not from the original start time for runs with continuation.
 294      * This has the advantage that one can find if conservation issues are
 295      * from the current run with the current settings on the current hardware.
 296      */
 297     void printEnergyConservation(FILE* fplog, int simulationPart, bool usingMdIntegrator) const;
 298
 299 private:
 300     //! Timestep
 301     double delta_t_ = 0;
 302
 303     //! Structure to store energy components and their running averages
 304     t_ebin* ebin_ = nullptr;
 305
 306     //! Is the periodic box triclinic
 307     bool bTricl_ = false;
 308     //! NHC trotter is used
 309     bool bNHC_trotter_ = false;
 310     //! If Nose-Hoover chains should be printed
 311     bool bPrintNHChains_ = false;
 312     //! If MTTK integrator was used
 313     bool bMTTK_ = false;
 314
 315     //! Temperature control scheme
 316     TemperatureCoupling etc_ = TemperatureCoupling::No;
 317
 318     //! Which of the main energy terms should be printed
 319     bool bEner_[F_NRE] = { false };
 320     //! Index for main energy terms
 321     int ie_ = 0;
 322     //! Number of energy terms from F_NRE list to be saved (i.e. number of 'true' in bEner)
 323     int f_nre_ = 0;
 324
 325     //! Index for constraints RMSD
 326     int iconrmsd_ = 0;
 327     /* !\brief How many constraints RMSD terms there are.
 328      * Usually 1 if LINCS is used and 0 otherwise)
 329      * nCrmsd > 0 indicates when constraints RMSD is saves, hence no boolean
 330      */
 331     int nCrmsd_ = 0;
 332
 333     //! Is the periodic box dynamic
 334     bool bDynBox_ = false;
 335     //! Index for box dimensions
 336     int ib_ = 0;
 337     //! Index for box volume
 338     int ivol_ = 0;
 339     //! Index for density
 340     int idens_ = 0;
 341     //! Triclinic box and not a rerun
 342     bool bDiagPres_ = false;
 343     //! Reference pressure, averaged over dimensions
 344     real ref_p_ = 0.0;
 345     //! Index for thermodynamic work (pV)
 346     int ipv_ = 0;
 347     //! Index for entalpy (pV + total energy)
 348     int ienthalpy_ = 0;
 349
 350     /*! \brief If the constraints virial should be printed.
 351      * Can only be true if "GMX_CONSTRAINTVIR" environmental variable is set */
 352     bool bConstrVir_ = false;
 353     //! Index for constrains virial
 354     int isvir_ = 0;
 355     //! Index for force virial
 356     int ifvir_ = 0;
 357
 358     //! If we have pressure computed
 359     bool bPres_ = false;
 360     //! Index for total virial
 361     int ivir_ = 0;
 362     //! Index for pressure
 363     int ipres_ = 0;
 364     /*! \brief Index for surface tension
 365      * [(pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ]]*/
 366     int isurft_ = 0;
 367
 368     //! Pressure control scheme
 369     PressureCoupling epc_ = PressureCoupling::No;
 370     //! Index for velocity of the box borders
 371     int ipc_ = 0;
 372
 373     //! If dipole was calculated
 374     bool bMu_ = false;
 375     //! Index for the dipole
 376     int imu_ = 0;
 377
 378     //! Index for coseine acceleration used for viscocity calculation
 379     int ivcos_ = 0;
 380     //! Index for viscocity
 381     int ivisc_ = 0;
 382
 383     //! Which energy terms from NonBondedEnergyTerms list should be printed in group-to-group block
 384     gmx::EnumerationArray<NonBondedEnergyTerms, bool> bEInd_;
 385     //! Number of energy terms to be printed (these, for which bEInd[] == true)
 386     int nEc_ = 0;
 387     //! Number of energy output groups
 388     int nEg_ = 0;
 389     //! Number of intergroup energy sets to be printed for each energy term (nE = (nEg*(nEg+1))/2)
 390     int nE_ = 0;
 391     //! Indexes for integroup energy sets (each set with nEc energies)
 392     std::vector<int> igrp_;
 393
 394     //! Number of temperature coupling groups
 395     int nTC_ = 0;
 396     //! Index for temperature
 397     int itemp_ = 0;
 398
 399     //! Number of Nose-Hoover chains
 400     int nNHC_ = 0;
 401     //! Number of temperature coupling coefficients in case of NH Chains
 402     int mde_n_ = 0;
 403     //! Index for temperature coupling coefficient in case of NH chains
 404     int itc_ = 0;
 405
 406     //! Number of temperature coupling terms if the temperature control is dealt by barostat (MTTK case)
 407     int nTCP_ = 0;
 408     //! Scalling factors for temperaturs control in MTTK
 409     int mdeb_n_ = 0;
 410     //! Index for scalling factor of MTTK
 411     int itcb_ = 0;
 412
 413     //! Array to accumulate values during update
 414     std::vector<real> tmp_r_;
 415
 416     //! The dhdl.xvg output file
 417     FILE* fp_dhdl_ = nullptr;
 418     //! Energy components for dhdl.xvg output
 419     std::vector<double> dE_;
 420     //! The delta U components (raw data + histogram)
 421     std::unique_ptr<t_mde_delta_h_coll> dhc_;
 422     //! Temperatures for simulated tempering groups
 423     std::vector<real> temperatures_;
 424
 425     //! For tracking the conserved or total energy
 426     std::unique_ptr<EnergyDriftTracker> conservedEnergyTracker_;
 427 };
 428
 429 } // namespace gmx
 430
 431 //! Open the dhdl file for output
 432 FILE* open_dhdl(const char* filename, const t_inputrec* ir, const gmx_output_env_t* oenv);
 433
 434 #endif