2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
35 /*! \libinternal \file
36 * \brief Header for the code that writes energy-like quantities.
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>
42 * \ingroup module_mdlib
45 #ifndef GMX_MDLIB_ENERGYOUTPUT_H
46 #define GMX_MDLIB_ENERGYOUTPUT_H
50 #include "gromacs/mdtypes/enerdata.h"
52 class energyhistory_t;
54 struct gmx_ekindata_t;
55 struct gmx_enerdata_t;
56 struct SimulationGroups;
58 struct gmx_output_env_t;
68 struct t_mde_delta_h_coll;
74 struct MdModulesNotifier;
77 //! \brief Printed names for intergroup energies
78 extern const char *egrp_nm[egNR+1];
80 /* \brief delta_h block type enum: the kinds of energies written out. */
83 //! Delta H BAR energy difference
85 //! dH/dlambda derivative
91 //! Expanded ensemble statistics
93 //! Total number of energy types in this enum
100 /* Energy output class
102 * This is the collection of energy averages collected during mdrun, and to
103 * be written out to the .edr file.
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)
113 /*! \brief Initiate MD energy bin
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.
123 EnergyOutput(ener_file * fp_ene,
124 const gmx_mtop_t * mtop,
125 const t_inputrec * ir,
126 const pull_t * pull_work,
129 const MdModulesNotifier &mdModulesNotifier);
133 /*! \brief Update the averaging structures.
135 * Called every step on which the thermodynamic values are evaluated.
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).
154 void addDataAtEnergyStep(bool bDoDHDL,
158 const gmx_enerdata_t *enerd,
159 const t_state *state,
161 const t_expanded *expand,
162 const matrix lastbox,
167 const gmx_ekindata_t *ekind,
169 const gmx::Constraints *constr);
171 /*! \brief Update the data averaging structure counts.
173 * Updates the number of steps, the values have not being computed.
175 void recordNonEnergyStep();
177 /*! \brief Writes current quantites to log and energy files.
179 * Prints current values of energies, pressure, temperature, restraint
180 * data, etc. to energy output file and to the log file (if not nullptr).
182 * This function only does something useful when bEne || bDR || bOR || log.
184 * \todo Perhaps this responsibility should involve some other
185 * object visiting all the contributing objects.
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.
198 void printStepToEnergyFile(ener_file *fp_ene, bool bEne, bool bDR, bool bOR,
200 int64_t step, double time,
204 /*! \brief Print reference temperatures for annealing groups.
206 * Nothing is done if log is nullptr.
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).
213 void printAnnealingTemperatures(FILE *log, SimulationGroups *groups, t_grpopts *opts);
215 /*! \brief Prints average values to log file.
217 * This is called at the end of the simulation run to print accumulated average values.
218 * Nothing it done if log is nullptr.
220 * \param[in] log Where to print.
221 * \param[in] groups Atom groups.
223 void printAverages(FILE *log, const SimulationGroups *groups);
225 /*! \brief Get the number of thermodynamic terms recorded.
227 * \todo Refactor this to return the expected output size,
228 * rather than exposing the implementation details about
229 * thermodynamic terms.
231 * \returns Number of thermodynamic terms.
233 int numEnergyTerms() const;
235 /*! \brief Fill the energyhistory_t data.
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.
242 * \param[out] enerhist Energy history data structure.
244 void fillEnergyHistory(energyhistory_t * enerhist) const;
246 /*! \brief Restore from energyhistory_t data.
248 * \param[in] enerhist Energy history data structure.
250 void restoreFromEnergyHistory(const energyhistory_t &enerhist);
252 //! Print an output header to the log file.
253 void printHeader(FILE *log, int64_t steps, double time);
259 //! Structure to store energy components and their running averages
260 t_ebin *ebin_ = nullptr;
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
271 //! Temperature control scheme
274 //! Which of the main energy terms should be printed
275 bool bEner_[F_NRE] = {false};
276 //! Index for main energy terms
278 //! Number of energy terms from F_NRE list to be saved (i.e. number of 'true' in bEner)
281 //! Index for constraints RMSD
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
289 //! Is the periodic box dynamic
290 bool bDynBox_ = false;
291 //! Index for box dimensions
293 //! Index for box volume
295 //! Index for density
297 //! Triclinic box and not a rerun
298 bool bDiagPres_ = false;
299 //! Reference pressure, averaged over dimensions
301 //! Index for thermodynamic work (pV)
303 //! Index for entalpy (pV + total energy)
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
311 //! Index for force virial
314 //! If we have pressure computed
316 //! Index for total virial
318 //! Index for pressure
320 /*! \brief Index for surface tension
321 * [(pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ]]*/
324 //! Pressure control scheme
326 //! Index for velocity of the box borders
329 //! If dipole was calculated
331 //! Index for the dipole
334 //! Index for coseine acceleration used for viscocity calculation
336 //! Index for viscocity
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)
343 //! Number of energy output groups
345 //! Number of intergroup energy sets to be printed for each energy term (nE = (nEg*(nEg+1))/2)
347 //! Indexes for integroup energy sets (each set with nEc energies)
348 int *igrp_ = nullptr;
350 //! Number of temperature coupling groups
352 //! Index for temperature
355 //! Number of Nose-Hoover chains
357 //! Number of temperature coupling coefficients in case of NH Chains
359 //! Index for temperature coupling coefficient in case of NH chains
362 //! Number of temperature coupling terms if the temperature control is dealt by barostat (MTTK case)
364 //! Scalling factors for temperaturs control in MTTK
366 //! Index for scalling factor of MTTK
369 //! Number of acceleration groups
371 //! Index for group velocities
374 //! Array to accumulate values during update
375 real *tmp_r_ = nullptr;
376 //! Array to accumulate values during update
377 rvec *tmp_v_ = nullptr;
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;
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);