src/gromacs/mdlib/energyoutput.h

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