2 * This file is part of the GROMACS molecular simulation package.
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.
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
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"
59 class energyhistory_t;
62 struct gmx_enerdata_t;
63 struct SimulationGroups;
65 struct gmx_output_env_t;
74 struct t_mde_delta_h_coll;
80 struct MDModulesNotifiers;
81 enum class StartingBehavior;
84 //! \brief Printed names for intergroup energies
85 const char* enumValueToString(NonBondedEnergyTerms enumValue);
87 /* \brief delta_h block type enum: the kinds of energies written out. */
90 //! Delta H BAR energy difference
92 //! dH/dlambda derivative
98 //! Expanded ensemble statistics
100 //! Total number of energy types in this enum
106 class EnergyDriftTracker;
109 * \brief Arrays connected to Pressure and Temperature coupling
111 struct PTCouplingArrays
113 //! Box velocities for Parrinello-Rahman P-coupling.
115 //! Nose-Hoover coordinates.
116 ArrayRef<const double> nosehoover_xi;
117 //! Nose-Hoover velocities.
118 ArrayRef<const double> nosehoover_vxi;
119 //! Pressure Nose-Hoover coordinates.
120 ArrayRef<const double> nhpres_xi;
121 //! Pressure Nose-Hoover velocities.
122 ArrayRef<const double> nhpres_vxi;
125 /* Energy output class
127 * This is the collection of energy averages collected during mdrun, and to
128 * be written out to the .edr file.
130 * \todo Use more std containers.
131 * \todo Write free-energy output also to energy file (after adding more tests)
136 /*! \brief Initiate MD energy bin
138 * \param[in] fp_ene Energy output file.
139 * \param[in] mtop Topology.
140 * \param[in] inputrec Input parameters.
141 * \param[in] pull_work Pulling simulations data
142 * \param[in] fp_dhdl FEP file.
143 * \param[in] isRerun Is this is a rerun instead of the simulations.
144 * \param[in] startingBehavior Run starting behavior.
145 * \param[in] simulationsShareState Tells whether the physical state is shared over simulations
146 * \param[in] mdModulesNotifiers Notifications to MD modules.
148 EnergyOutput(ener_file* fp_ene,
149 const gmx_mtop_t& mtop,
150 const t_inputrec& inputrec,
151 const pull_t* pull_work,
154 StartingBehavior startingBehavior,
155 bool simulationsShareState,
156 const MDModulesNotifiers& mdModulesNotifiers);
160 /*! \brief Update the averaging structures.
162 * Called every step on which the thermodynamic values are evaluated.
164 * \param[in] bDoDHDL Whether the FEP is enabled.
165 * \param[in] bSum If this stepshould be recorded to compute sums and averages.
166 * \param[in] time Current simulation time.
167 * \param[in] tmass Total mass
168 * \param[in] enerd Energy data object.
169 * \param[in] fep FEP data.
170 * \param[in] lastbox PBC data.
171 * \param[in] ptCouplingArrays Arrays connected to pressure and temperature coupling.
172 * \param[in] fep_state The current alchemical state we are in.
173 * \param[in] vir Total virial.
174 * \param[in] pres Pressure.
175 * \param[in] ekind Kinetic energy data.
176 * \param[in] mu_tot Total dipole.
177 * \param[in] constr Constraints object to get RMSD from (for LINCS only).
179 void addDataAtEnergyStep(bool bDoDHDL,
183 const gmx_enerdata_t* enerd,
185 const matrix lastbox,
186 PTCouplingArrays ptCouplingArrays,
190 const gmx_ekindata_t* ekind,
192 const gmx::Constraints* constr);
194 /*! \brief Update the data averaging structure counts.
196 * Updates the number of steps, the values have not being computed.
198 void recordNonEnergyStep();
200 /*! \brief Writes current quantites to log and energy files.
202 * Prints current values of energies, pressure, temperature, restraint
203 * data, etc. to energy output file and to the log file (if not nullptr).
205 * This function only does something useful when bEne || bDR || bOR || log.
207 * \todo Perhaps this responsibility should involve some other
208 * object visiting all the contributing objects.
210 * \param[in] fp_ene Energy file for the output.
211 * \param[in] bEne If it is a step for energy output or last step.
212 * \param[in] bDR If it is a step of writing distance restraints.
213 * \param[in] bOR If it is a step of writing orientation restraints.
214 * \param[in] log Pointer to the log file.
215 * \param[in] step Current step.
216 * \param[in] time Current simulation time.
217 * \param[in] fcd Bonded force computation data,
218 * including orientation and distance restraints.
219 * \param[in] awh AWH data.
221 void printStepToEnergyFile(ener_file* fp_ene,
231 /*! \brief Print reference temperatures for annealing groups.
233 * Nothing is done if log is nullptr.
235 * \param[in] log Log file to print to.
236 * \param[in] groups Information on atom groups.
237 * \param[in] opts Atom temperature coupling groups options
238 * (annealing is done by groups).
240 static void printAnnealingTemperatures(FILE* log, const SimulationGroups* groups, const t_grpopts* opts);
242 /*! \brief Prints average values to log file.
244 * This is called at the end of the simulation run to print accumulated average values.
245 * Nothing it done if log is nullptr.
247 * \param[in] log Where to print.
248 * \param[in] groups Atom groups.
250 void printAverages(FILE* log, const SimulationGroups* groups);
252 /*! \brief Get the number of thermodynamic terms recorded.
254 * \todo Refactor this to return the expected output size,
255 * rather than exposing the implementation details about
256 * thermodynamic terms.
258 * \returns Number of thermodynamic terms.
260 int numEnergyTerms() const;
262 /*! \brief Fill the energyhistory_t data.
264 * Between .edr writes, the averages are history dependent,
265 * and that history needs to be retained in checkpoints.
266 * These functions set/read the energyhistory_t class
267 * that is written to checkpoints.
269 * \param[out] enerhist Energy history data structure.
271 void fillEnergyHistory(energyhistory_t* enerhist) const;
273 /*! \brief Restore from energyhistory_t data.
275 * \param[in] enerhist Energy history data structure.
277 void restoreFromEnergyHistory(const energyhistory_t& enerhist);
279 //! Print an output header to the log file.
280 static void printHeader(FILE* log, int64_t steps, double time);
282 /*! \brief Print conserved energy drift message to \p fplog
284 * Note that this is only over the current run (times are printed),
285 * this is not from the original start time for runs with continuation.
286 * This has the advantage that one can find if conservation issues are
287 * from the current run with the current settings on the current hardware.
289 void printEnergyConservation(FILE* fplog, int simulationPart, bool usingMdIntegrator) const;
295 //! Structure to store energy components and their running averages
296 t_ebin* ebin_ = nullptr;
298 //! Is the periodic box triclinic
299 bool bTricl_ = false;
300 //! NHC trotter is used
301 bool bNHC_trotter_ = false;
302 //! If Nose-Hoover chains should be printed
303 bool bPrintNHChains_ = false;
304 //! If MTTK integrator was used
307 //! Temperature control scheme
308 TemperatureCoupling etc_ = TemperatureCoupling::No;
310 //! Which of the main energy terms should be printed
311 bool bEner_[F_NRE] = { false };
312 //! Index for main energy terms
314 //! Number of energy terms from F_NRE list to be saved (i.e. number of 'true' in bEner)
317 //! Index for constraints RMSD
319 /* !\brief How many constraints RMSD terms there are.
320 * Usually 1 if LINCS is used and 0 otherwise)
321 * nCrmsd > 0 indicates when constraints RMSD is saves, hence no boolean
325 //! Is the periodic box dynamic
326 bool bDynBox_ = false;
327 //! Index for box dimensions
329 //! Index for box volume
331 //! Index for density
333 //! Triclinic box and not a rerun
334 bool bDiagPres_ = false;
335 //! Reference pressure, averaged over dimensions
337 //! Index for thermodynamic work (pV)
339 //! Index for entalpy (pV + total energy)
342 //! If we have pressure computed
344 //! Index for total virial
346 //! Index for pressure
348 /*! \brief Index for surface tension
349 * [(pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ]]*/
352 //! Pressure control scheme
353 PressureCoupling epc_ = PressureCoupling::No;
354 //! Index for velocity of the box borders
357 //! If dipole was calculated
359 //! Index for the dipole
362 //! Index for coseine acceleration used for viscocity calculation
364 //! Index for viscocity
367 //! Which energy terms from NonBondedEnergyTerms list should be printed in group-to-group block
368 gmx::EnumerationArray<NonBondedEnergyTerms, bool> bEInd_;
369 //! Number of energy terms to be printed (these, for which bEInd[] == true)
371 //! Number of energy output groups
373 //! Number of intergroup energy sets to be printed for each energy term (nE = (nEg*(nEg+1))/2)
375 //! Indexes for integroup energy sets (each set with nEc energies)
376 std::vector<int> igrp_;
378 //! Number of temperature coupling groups
380 //! Index for temperature
383 //! Number of Nose-Hoover chains
385 //! Number of temperature coupling coefficients in case of NH Chains
387 //! Index for temperature coupling coefficient in case of NH chains
390 //! Number of temperature coupling terms if the temperature control is dealt by barostat (MTTK case)
392 //! Scalling factors for temperaturs control in MTTK
394 //! Index for scalling factor of MTTK
397 //! Array to accumulate values during update
398 std::vector<real> tmp_r_;
400 //! The dhdl.xvg output file
401 FILE* fp_dhdl_ = nullptr;
402 //! Whether the free-energy lambda moves dynamically between lambda states
403 bool haveFepLambdaMoves_;
404 //! Energy components for dhdl.xvg output
405 std::vector<double> dE_;
406 //! The delta U components (raw data + histogram)
407 std::unique_ptr<t_mde_delta_h_coll> dhc_;
408 //! Temperatures for simulated tempering groups
409 std::vector<real> temperatures_;
411 //! For tracking the conserved or total energy
412 std::unique_ptr<EnergyDriftTracker> conservedEnergyTracker_;
417 //! Open the dhdl file for output
418 FILE* open_dhdl(const char* filename, const t_inputrec* ir, const gmx_output_env_t* oenv);