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