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