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