src/gromacs/mdlib/energyoutput.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  * \brief Defines code that writes energy-like quantities.
  39  *
  40  * \author Mark Abraham <mark.j.abraham@gmail.com>
  41  * \author Paul Bauer <paul.bauer.q@gmail.com>
  42  * \author Artem Zhmurov <zhmurov@gmail.com>
  43  *
  44  * \ingroup module_mdlib
  45  */
  46 #include "gmxpre.h"
  47
  48 #include "energyoutput.h"
  49
  50 #include <cfloat>
  51 #include <cstdlib>
  52 #include <cstring>
  53
  54 #include <array>
  55 #include <string>
  56
  57 #include "gromacs/awh/awh.h"
  58 #include "gromacs/fileio/enxio.h"
  59 #include "gromacs/fileio/gmxfio.h"
  60 #include "gromacs/fileio/xvgr.h"
  61 #include "gromacs/gmxlib/network.h"
  62 #include "gromacs/listed_forces/disre.h"
  63 #include "gromacs/listed_forces/orires.h"
  64 #include "gromacs/math/functions.h"
  65 #include "gromacs/math/units.h"
  66 #include "gromacs/math/vec.h"
  67 #include "gromacs/mdlib/constr.h"
  68 #include "gromacs/mdlib/ebin.h"
  69 #include "gromacs/mdlib/mdebin_bar.h"
  70 #include "gromacs/mdtypes/energyhistory.h"
  71 #include "gromacs/mdtypes/fcdata.h"
  72 #include "gromacs/mdtypes/group.h"
  73 #include "gromacs/mdtypes/inputrec.h"
  74 #include "gromacs/mdtypes/md_enums.h"
  75 #include "gromacs/mdtypes/state.h"
  76 #include "gromacs/pbcutil/pbc.h"
  77 #include "gromacs/pulling/pull.h"
  78 #include "gromacs/topology/mtop_util.h"
  79 #include "gromacs/trajectory/energyframe.h"
  80 #include "gromacs/utility/arraysize.h"
  81 #include "gromacs/utility/fatalerror.h"
  82 #include "gromacs/utility/smalloc.h"
  83 #include "gromacs/utility/stringutil.h"
  84
  85 //! Labels for energy file quantities
  86 //! \{
  87 static const char                 *conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
  88
  89 static std::array<const char *, 3> boxs_nm = { "Box-X", "Box-Y", "Box-Z" };
  90
  91 static std::array<const char *, 6> tricl_boxs_nm = {
  92     "Box-XX", "Box-YY", "Box-ZZ",
  93     "Box-YX", "Box-ZX", "Box-ZY"
  94 };
  95
  96 static const char                 *vol_nm[] = { "Volume" };
  97
  98 static const char                 *dens_nm[] = {"Density" };
  99
 100 static const char                 *pv_nm[] = {"pV" };
 101
 102 static const char                 *enthalpy_nm[] = {"Enthalpy" };
 103
 104 static std::array<const char *, 6> boxvel_nm = {
 105     "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
 106     "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
 107 };
 108
 109 #define NBOXS asize(boxs_nm)
 110 #define NTRICLBOXS asize(tricl_boxs_nm)
 111
 112 const char *egrp_nm[egNR+1] = {
 113     "Coul-SR", "LJ-SR", "Buck-SR",
 114     "Coul-14", "LJ-14", nullptr
 115 };
 116 //! \}
 117
 118 namespace gmx
 119 {
 120
 121 /*! \brief 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 EnergyOutput::EnergyOutput(ener_file        *fp_ene,
 131                            const gmx_mtop_t *mtop,
 132                            const t_inputrec *ir,
 133                            const pull_t     *pull_work,
 134                            FILE             *fp_dhdl,
 135                            bool              isRerun)
 136 {
 137     const char                *ener_nm[F_NRE];
 138     static const char         *vir_nm[] = {
 139         "Vir-XX", "Vir-XY", "Vir-XZ",
 140         "Vir-YX", "Vir-YY", "Vir-YZ",
 141         "Vir-ZX", "Vir-ZY", "Vir-ZZ"
 142     };
 143     static const char         *sv_nm[] = {
 144         "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
 145         "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
 146         "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
 147     };
 148     static const char         *fv_nm[] = {
 149         "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
 150         "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
 151         "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
 152     };
 153     static const char         *pres_nm[] = {
 154         "Pres-XX", "Pres-XY", "Pres-XZ",
 155         "Pres-YX", "Pres-YY", "Pres-YZ",
 156         "Pres-ZX", "Pres-ZY", "Pres-ZZ"
 157     };
 158     static const char         *surft_nm[] = {
 159         "#Surf*SurfTen"
 160     };
 161     static const char         *mu_nm[] = {
 162         "Mu-X", "Mu-Y", "Mu-Z"
 163     };
 164     static const char         *vcos_nm[] = {
 165         "2CosZ*Vel-X"
 166     };
 167     static const char         *visc_nm[] = {
 168         "1/Viscosity"
 169     };
 170     static const char         *baro_nm[] = {
 171         "Barostat"
 172     };
 173
 174     const SimulationGroups    *groups;
 175     char                     **gnm;
 176     char                       buf[256];
 177     const char                *bufi;
 178     int                        i, j, ni, nj, n, k, kk, ncon, nset;
 179     bool                       bBHAM, b14;
 180
 181     if (EI_DYNAMICS(ir->eI))
 182     {
 183         delta_t_ = ir->delta_t;
 184     }
 185     else
 186     {
 187         delta_t_ = 0;
 188     }
 189
 190     groups = &mtop->groups;
 191
 192     bBHAM = (mtop->ffparams.numTypes() > 0) && (mtop->ffparams.functype[0] == F_BHAM);
 193     b14   = (gmx_mtop_ftype_count(mtop, F_LJ14) > 0 ||
 194              gmx_mtop_ftype_count(mtop, F_LJC14_Q) > 0);
 195
 196     ncon           = gmx_mtop_ftype_count(mtop, F_CONSTR);
 197     nset           = gmx_mtop_ftype_count(mtop, F_SETTLE);
 198     bool bConstr   = (ncon > 0 || nset > 0) && !isRerun;
 199     bConstrVir_ = false;
 200     nCrmsd_     = 0;
 201     if (bConstr)
 202     {
 203         if (ncon > 0 && ir->eConstrAlg == econtLINCS)
 204         {
 205             nCrmsd_ = 1;
 206         }
 207         bConstrVir_ = (getenv("GMX_CONSTRAINTVIR") != nullptr);
 208     }
 209     else
 210     {
 211         nCrmsd_ = 0;
 212     }
 213
 214     /* Energy monitoring */
 215     for (i = 0; i < egNR; i++)
 216     {
 217         bEInd_[i] = false;
 218     }
 219
 220     for (i = 0; i < F_NRE; i++)
 221     {
 222         bEner_[i] = false;
 223         if (isRerun &&
 224             (i == F_EKIN || i == F_ETOT || i == F_ECONSERVED ||
 225              i == F_TEMP || i == F_PDISPCORR || i == F_PRES))
 226         {
 227             continue;
 228         }
 229         if (i == F_LJ)
 230         {
 231             bEner_[i] = !bBHAM;
 232         }
 233         else if (i == F_BHAM)
 234         {
 235             bEner_[i] = bBHAM;
 236         }
 237         else if (i == F_EQM)
 238         {
 239             bEner_[i] = ir->bQMMM;
 240         }
 241         else if (i == F_RF_EXCL)
 242         {
 243             bEner_[i] = (EEL_RF(ir->coulombtype) && ir->cutoff_scheme == ecutsGROUP);
 244         }
 245         else if (i == F_COUL_RECIP)
 246         {
 247             bEner_[i] = EEL_FULL(ir->coulombtype);
 248         }
 249         else if (i == F_LJ_RECIP)
 250         {
 251             bEner_[i] = EVDW_PME(ir->vdwtype);
 252         }
 253         else if (i == F_LJ14)
 254         {
 255             bEner_[i] = b14;
 256         }
 257         else if (i == F_COUL14)
 258         {
 259             bEner_[i] = b14;
 260         }
 261         else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
 262         {
 263             bEner_[i] = false;
 264         }
 265         else if ((i == F_DVDL_COUL && ir->fepvals->separate_dvdl[efptCOUL]) ||
 266                  (i == F_DVDL_VDW  && ir->fepvals->separate_dvdl[efptVDW]) ||
 267                  (i == F_DVDL_BONDED && ir->fepvals->separate_dvdl[efptBONDED]) ||
 268                  (i == F_DVDL_RESTRAINT && ir->fepvals->separate_dvdl[efptRESTRAINT]) ||
 269                  (i == F_DKDL && ir->fepvals->separate_dvdl[efptMASS]) ||
 270                  (i == F_DVDL && ir->fepvals->separate_dvdl[efptFEP]))
 271         {
 272             bEner_[i] = (ir->efep != efepNO);
 273         }
 274         else if ((interaction_function[i].flags & IF_VSITE) ||
 275                  (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
 276         {
 277             bEner_[i] = false;
 278         }
 279         else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES)  || (i == F_EQM))
 280         {
 281             bEner_[i] = true;
 282         }
 283         else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
 284         {
 285             bEner_[i] = EI_DYNAMICS(ir->eI);
 286         }
 287         else if (i == F_DISPCORR || i == F_PDISPCORR)
 288         {
 289             bEner_[i] = (ir->eDispCorr != edispcNO);
 290         }
 291         else if (i == F_DISRESVIOL)
 292         {
 293             bEner_[i] = (gmx_mtop_ftype_count(mtop, F_DISRES) > 0);
 294         }
 295         else if (i == F_ORIRESDEV)
 296         {
 297             bEner_[i] = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
 298         }
 299         else if (i == F_CONNBONDS)
 300         {
 301             bEner_[i] = false;
 302         }
 303         else if (i == F_COM_PULL)
 304         {
 305             bEner_[i] = ((ir->bPull && pull_have_potential(pull_work)) ||
 306                          ir->bRot);
 307         }
 308         else if (i == F_ECONSERVED)
 309         {
 310             bEner_[i] = (integratorHasConservedEnergyQuantity(ir));
 311         }
 312         else
 313         {
 314             bEner_[i] = (gmx_mtop_ftype_count(mtop, i) > 0);
 315         }
 316     }
 317
 318     f_nre_ = 0;
 319     for (i = 0; i < F_NRE; i++)
 320     {
 321         if (bEner_[i])
 322         {
 323             ener_nm[f_nre_] = interaction_function[i].longname;
 324             f_nre_++;
 325         }
 326     }
 327
 328     epc_            = isRerun ? epcNO : ir->epc;
 329     bDiagPres_      = !TRICLINIC(ir->ref_p) && !isRerun;
 330     ref_p_          = (ir->ref_p[XX][XX]+ir->ref_p[YY][YY]+ir->ref_p[ZZ][ZZ])/DIM;
 331     bTricl_         = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
 332     bDynBox_        = inputrecDynamicBox(ir);
 333     etc_            = isRerun ? etcNO : ir->etc;
 334     bNHC_trotter_   = inputrecNvtTrotter(ir) && !isRerun;
 335     bPrintNHChains_ = ir->bPrintNHChains && !isRerun;
 336     bMTTK_          = (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)) && !isRerun;
 337     bMu_            = inputrecNeedMutot(ir);
 338     bPres_          = !isRerun;
 339
 340     ebin_  = mk_ebin();
 341     /* Pass NULL for unit to let get_ebin_space determine the units
 342      * for interaction_function[i].longname
 343      */
 344     ie_    = get_ebin_space(ebin_, f_nre_, ener_nm, nullptr);
 345     if (nCrmsd_)
 346     {
 347         /* This should be called directly after the call for ie_,
 348          * such that iconrmsd_ follows directly in the list.
 349          */
 350         iconrmsd_ = get_ebin_space(ebin_, nCrmsd_, conrmsd_nm, "");
 351     }
 352     if (bDynBox_)
 353     {
 354         ib_    = get_ebin_space(ebin_,
 355                                 bTricl_ ? tricl_boxs_nm.size() : boxs_nm.size(),
 356                                 bTricl_ ? tricl_boxs_nm.data() : boxs_nm.data(),
 357                                 unit_length);
 358         ivol_  = get_ebin_space(ebin_, 1, vol_nm,  unit_volume);
 359         idens_ = get_ebin_space(ebin_, 1, dens_nm, unit_density_SI);
 360         if (bDiagPres_)
 361         {
 362             ipv_       = get_ebin_space(ebin_, 1, pv_nm,   unit_energy);
 363             ienthalpy_ = get_ebin_space(ebin_, 1, enthalpy_nm,   unit_energy);
 364         }
 365     }
 366     if (bConstrVir_)
 367     {
 368         isvir_ = get_ebin_space(ebin_, asize(sv_nm), sv_nm, unit_energy);
 369         ifvir_ = get_ebin_space(ebin_, asize(fv_nm), fv_nm, unit_energy);
 370     }
 371     if (bPres_)
 372     {
 373         ivir_   = get_ebin_space(ebin_, asize(vir_nm), vir_nm, unit_energy);
 374         ipres_  = get_ebin_space(ebin_, asize(pres_nm), pres_nm, unit_pres_bar);
 375         isurft_ = get_ebin_space(ebin_, asize(surft_nm), surft_nm,
 376                                  unit_surft_bar);
 377     }
 378     if (epc_ == epcPARRINELLORAHMAN || epc_ == epcMTTK)
 379     {
 380         ipc_ = get_ebin_space(ebin_, bTricl_ ? boxvel_nm.size() : DIM,
 381                               boxvel_nm.data(), unit_vel);
 382     }
 383     if (bMu_)
 384     {
 385         imu_    = get_ebin_space(ebin_, asize(mu_nm), mu_nm, unit_dipole_D);
 386     }
 387     if (ir->cos_accel != 0)
 388     {
 389         ivcos_ = get_ebin_space(ebin_, asize(vcos_nm), vcos_nm, unit_vel);
 390         ivisc_ = get_ebin_space(ebin_, asize(visc_nm), visc_nm,
 391                                 unit_invvisc_SI);
 392     }
 393
 394     /* Energy monitoring */
 395     for (i = 0; i < egNR; i++)
 396     {
 397         bEInd_[i] = false;
 398     }
 399     bEInd_[egCOULSR] = true;
 400     bEInd_[egLJSR  ] = true;
 401
 402     if (bBHAM)
 403     {
 404         bEInd_[egLJSR]   = false;
 405         bEInd_[egBHAMSR] = true;
 406     }
 407     if (b14)
 408     {
 409         bEInd_[egLJ14]   = true;
 410         bEInd_[egCOUL14] = true;
 411     }
 412     nEc_ = 0;
 413     for (i = 0; (i < egNR); i++)
 414     {
 415         if (bEInd_[i])
 416         {
 417             nEc_++;
 418         }
 419     }
 420     n         = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
 421     nEg_      = n;
 422     nE_       = (n*(n+1))/2;
 423
 424     snew(igrp_, nE_);
 425     if (nE_ > 1)
 426     {
 427         n = 0;
 428         snew(gnm, nEc_);
 429         for (k = 0; (k < nEc_); k++)
 430         {
 431             snew(gnm[k], STRLEN);
 432         }
 433         for (i = 0; (i < gmx::ssize(groups->groups[SimulationAtomGroupType::EnergyOutput])); i++)
 434         {
 435             ni = groups->groups[SimulationAtomGroupType::EnergyOutput][i];
 436             for (j = i; (j < gmx::ssize(groups->groups[SimulationAtomGroupType::EnergyOutput])); j++)
 437             {
 438                 nj = groups->groups[SimulationAtomGroupType::EnergyOutput][j];
 439                 for (k = kk = 0; (k < egNR); k++)
 440                 {
 441                     if (bEInd_[k])
 442                     {
 443                         sprintf(gnm[kk], "%s:%s-%s", egrp_nm[k],
 444                                 *(groups->groupNames[ni]), *(groups->groupNames[nj]));
 445                         kk++;
 446                     }
 447                 }
 448                 igrp_[n] = get_ebin_space(ebin_, nEc_,
 449                                           gnm, unit_energy);
 450                 n++;
 451             }
 452         }
 453         for (k = 0; (k < nEc_); k++)
 454         {
 455             sfree(gnm[k]);
 456         }
 457         sfree(gnm);
 458
 459         if (n != nE_)
 460         {
 461             gmx_incons("Number of energy terms wrong");
 462         }
 463     }
 464
 465     nTC_  = isRerun ? 0 : groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
 466     nNHC_ = ir->opts.nhchainlength; /* shorthand for number of NH chains */
 467     if (bMTTK_)
 468     {
 469         nTCP_ = 1;  /* assume only one possible coupling system for barostat
 470                             for now */
 471     }
 472     else
 473     {
 474         nTCP_ = 0;
 475     }
 476     if (etc_ == etcNOSEHOOVER)
 477     {
 478         if (bNHC_trotter_)
 479         {
 480             mde_n_ = 2*nNHC_*nTC_;
 481         }
 482         else
 483         {
 484             mde_n_ = 2*nTC_;
 485         }
 486         if (epc_ == epcMTTK)
 487         {
 488             mdeb_n_ = 2*nNHC_*nTCP_;
 489         }
 490     }
 491     else
 492     {
 493         mde_n_  = nTC_;
 494         mdeb_n_ = 0;
 495     }
 496
 497     snew(tmp_r_, mde_n_);
 498     // TODO redo the group name memory management to make it more clear
 499     char **grpnms;
 500     snew(grpnms, std::max(mde_n_, mdeb_n_)); // Just in case mdeb_n_ > mde_n_
 501
 502     for (i = 0; (i < nTC_); i++)
 503     {
 504         ni = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
 505         sprintf(buf, "T-%s", *(groups->groupNames[ni]));
 506         grpnms[i] = gmx_strdup(buf);
 507     }
 508     itemp_ = get_ebin_space(ebin_, nTC_, grpnms,
 509                             unit_temp_K);
 510     for (i = 0; i < nTC_; i++)
 511     {
 512         sfree(grpnms[i]);
 513     }
 514
 515     int allocated = 0;
 516     if (etc_ == etcNOSEHOOVER)
 517     {
 518         if (bPrintNHChains_)
 519         {
 520             if (bNHC_trotter_)
 521             {
 522                 for (i = 0; (i < nTC_); i++)
 523                 {
 524                     ni   = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
 525                     bufi = *(groups->groupNames[ni]);
 526                     for (j = 0; (j < nNHC_); j++)
 527                     {
 528                         sprintf(buf, "Xi-%d-%s", j, bufi);
 529                         grpnms[2*(i*nNHC_+j)] = gmx_strdup(buf);
 530                         sprintf(buf, "vXi-%d-%s", j, bufi);
 531                         grpnms[2*(i*nNHC_+j)+1] = gmx_strdup(buf);
 532                     }
 533                 }
 534                 itc_ = get_ebin_space(ebin_, mde_n_,
 535                                       grpnms, unit_invtime);
 536                 allocated = mde_n_;
 537                 if (bMTTK_)
 538                 {
 539                     for (i = 0; (i < nTCP_); i++)
 540                     {
 541                         bufi = baro_nm[0];  /* All barostat DOF's together for now. */
 542                         for (j = 0; (j < nNHC_); j++)
 543                         {
 544                             sprintf(buf, "Xi-%d-%s", j, bufi);
 545                             grpnms[2*(i*nNHC_+j)] = gmx_strdup(buf);
 546                             sprintf(buf, "vXi-%d-%s", j, bufi);
 547                             grpnms[2*(i*nNHC_+j)+1] = gmx_strdup(buf);
 548                         }
 549                     }
 550                     itcb_ = get_ebin_space(ebin_, mdeb_n_,
 551                                            grpnms, unit_invtime);
 552                     allocated = mdeb_n_;
 553                 }
 554             }
 555             else
 556             {
 557                 for (i = 0; (i < nTC_); i++)
 558                 {
 559                     ni   = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
 560                     bufi = *(groups->groupNames[ni]);
 561                     sprintf(buf, "Xi-%s", bufi);
 562                     grpnms[2*i] = gmx_strdup(buf);
 563                     sprintf(buf, "vXi-%s", bufi);
 564                     grpnms[2*i+1] = gmx_strdup(buf);
 565                 }
 566                 itc_ = get_ebin_space(ebin_, mde_n_,
 567                                       grpnms, unit_invtime);
 568                 allocated = mde_n_;
 569             }
 570         }
 571     }
 572     else if (etc_ == etcBERENDSEN || etc_ == etcYES ||
 573              etc_ == etcVRESCALE)
 574     {
 575         for (i = 0; (i < nTC_); i++)
 576         {
 577             ni = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
 578             sprintf(buf, "Lamb-%s", *(groups->groupNames[ni]));
 579             grpnms[i] = gmx_strdup(buf);
 580         }
 581         itc_        = get_ebin_space(ebin_, mde_n_, grpnms, "");
 582         allocated   = mde_n_;
 583     }
 584
 585     for (i = 0; i < allocated; i++)
 586     {
 587         sfree(grpnms[i]);
 588     }
 589     sfree(grpnms);
 590
 591     nU_ = groups->groups[SimulationAtomGroupType::Acceleration].size();
 592     snew(tmp_v_, nU_);
 593     if (nU_ > 1)
 594     {
 595         snew(grpnms, 3*nU_);
 596         for (i = 0; (i < nU_); i++)
 597         {
 598             ni = groups->groups[SimulationAtomGroupType::Acceleration][i];
 599             sprintf(buf, "Ux-%s", *(groups->groupNames[ni]));
 600             grpnms[3*i+XX] = gmx_strdup(buf);
 601             sprintf(buf, "Uy-%s", *(groups->groupNames[ni]));
 602             grpnms[3*i+YY] = gmx_strdup(buf);
 603             sprintf(buf, "Uz-%s", *(groups->groupNames[ni]));
 604             grpnms[3*i+ZZ] = gmx_strdup(buf);
 605         }
 606         iu_ = get_ebin_space(ebin_, 3*nU_, grpnms, unit_vel);
 607         for (i = 0; i < 3*nU_; i++)
 608         {
 609             sfree(grpnms[i]);
 610         }
 611         sfree(grpnms);
 612     }
 613
 614     if (fp_ene)
 615     {
 616         do_enxnms(fp_ene, &ebin_->nener, &ebin_->enm);
 617     }
 618
 619     /* check whether we're going to write dh histograms */
 620     dhc_ = nullptr;
 621     if (ir->fepvals->separate_dhdl_file == esepdhdlfileNO)
 622     {
 623         /* Currently dh histograms are only written with dynamics */
 624         if (EI_DYNAMICS(ir->eI))
 625         {
 626             snew(dhc_, 1);
 627
 628             mde_delta_h_coll_init(dhc_, ir);
 629         }
 630         fp_dhdl_ = nullptr;
 631         snew(dE_, ir->fepvals->n_lambda);
 632     }
 633     else
 634     {
 635         fp_dhdl_ = fp_dhdl;
 636         snew(dE_, ir->fepvals->n_lambda);
 637     }
 638     if (ir->bSimTemp)
 639     {
 640         int i;
 641         snew(temperatures_, ir->fepvals->n_lambda);
 642         numTemperatures_ = ir->fepvals->n_lambda;
 643         for (i = 0; i < ir->fepvals->n_lambda; i++)
 644         {
 645             temperatures_[i] = ir->simtempvals->temperatures[i];
 646         }
 647     }
 648     else
 649     {
 650         numTemperatures_ = 0;
 651     }
 652
 653 }
 654
 655 EnergyOutput::~EnergyOutput()
 656 {
 657     sfree(igrp_);
 658     sfree(tmp_r_);
 659     sfree(tmp_v_);
 660     done_ebin(ebin_);
 661     done_mde_delta_h_coll(dhc_);
 662     sfree(dE_);
 663     if (numTemperatures_ > 0)
 664     {
 665         sfree(temperatures_);
 666     }
 667 }
 668
 669 } // namespace gmx
 670
 671 /*! \brief Print a lambda vector to a string
 672  *
 673  * \param[in] fep                The inputrec's FEP input data
 674  * \param[in] i                  The index of the lambda vector
 675  * \param[in] get_native_lambda  Whether to print the native lambda
 676  * \param[in] get_names          Whether to print the names rather than the values
 677  * \param[in,out] str            The pre-allocated string buffer to print to.
 678  */
 679 static void print_lambda_vector(t_lambda *fep, int i,
 680                                 bool get_native_lambda, bool get_names,
 681                                 char *str)
 682 {
 683     int    j, k = 0;
 684     int    Nsep = 0;
 685
 686     for (j = 0; j < efptNR; j++)
 687     {
 688         if (fep->separate_dvdl[j])
 689         {
 690             Nsep++;
 691         }
 692     }
 693     str[0] = 0; /* reset the string */
 694     if (Nsep > 1)
 695     {
 696         str += sprintf(str, "("); /* set the opening parenthesis*/
 697     }
 698     for (j = 0; j < efptNR; j++)
 699     {
 700         if (fep->separate_dvdl[j])
 701         {
 702             if (!get_names)
 703             {
 704                 if (get_native_lambda && fep->init_lambda >= 0)
 705                 {
 706                     str += sprintf(str, "%.4f", fep->init_lambda);
 707                 }
 708                 else
 709                 {
 710                     str += sprintf(str, "%.4f", fep->all_lambda[j][i]);
 711                 }
 712             }
 713             else
 714             {
 715                 str += sprintf(str, "%s", efpt_singular_names[j]);
 716             }
 717             /* print comma for the next item */
 718             if (k < Nsep-1)
 719             {
 720                 str += sprintf(str, ", ");
 721             }
 722             k++;
 723         }
 724     }
 725     if (Nsep > 1)
 726     {
 727         /* and add the closing parenthesis */
 728         sprintf(str, ")");
 729     }
 730 }
 731
 732 FILE *open_dhdl(const char *filename, const t_inputrec *ir,
 733                 const gmx_output_env_t *oenv)
 734 {
 735     FILE       *fp;
 736     const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda",
 737     *lambdastate     = "\\lambda state";
 738     int         i, nsets, nsets_de, nsetsbegin;
 739     int         n_lambda_terms = 0;
 740     t_lambda   *fep            = ir->fepvals; /* for simplicity */
 741     t_expanded *expand         = ir->expandedvals;
 742     char        lambda_vec_str[STRLEN], lambda_name_str[STRLEN];
 743
 744     int         nsets_dhdl = 0;
 745     int         s          = 0;
 746     int         nsetsextend;
 747     bool        write_pV = false;
 748
 749     /* count the number of different lambda terms */
 750     for (i = 0; i < efptNR; i++)
 751     {
 752         if (fep->separate_dvdl[i])
 753         {
 754             n_lambda_terms++;
 755         }
 756     }
 757
 758     std::string title, label_x, label_y;
 759     if (fep->n_lambda == 0)
 760     {
 761         title   = gmx::formatString("%s", dhdl);
 762         label_x = gmx::formatString("Time (ps)");
 763         label_y = gmx::formatString("%s (%s %s)",
 764                                     dhdl, unit_energy, "[\\lambda]\\S-1\\N");
 765     }
 766     else
 767     {
 768         title   = gmx::formatString("%s and %s", dhdl, deltag);
 769         label_x = gmx::formatString("Time (ps)");
 770         label_y = gmx::formatString("%s and %s (%s %s)",
 771                                     dhdl, deltag, unit_energy, "[\\8l\\4]\\S-1\\N");
 772     }
 773     fp = gmx_fio_fopen(filename, "w+");
 774     xvgr_header(fp, title.c_str(), label_x, label_y, exvggtXNY, oenv);
 775
 776     std::string buf;
 777     if (!(ir->bSimTemp))
 778     {
 779         buf = gmx::formatString("T = %g (K) ", ir->opts.ref_t[0]);
 780     }
 781     if ((ir->efep != efepSLOWGROWTH) && (ir->efep != efepEXPANDED))
 782     {
 783         if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))
 784         {
 785             /* compatibility output */
 786             buf += gmx::formatString("%s = %.4f", lambda, fep->init_lambda);
 787         }
 788         else
 789         {
 790             print_lambda_vector(fep, fep->init_fep_state, true, false,
 791                                 lambda_vec_str);
 792             print_lambda_vector(fep, fep->init_fep_state, true, true,
 793                                 lambda_name_str);
 794             buf += gmx::formatString("%s %d: %s = %s",
 795                                      lambdastate, fep->init_fep_state,
 796                                      lambda_name_str, lambda_vec_str);
 797         }
 798     }
 799     xvgr_subtitle(fp, buf.c_str(), oenv);
 800
 801
 802     nsets_dhdl = 0;
 803     if (fep->dhdl_derivatives == edhdlderivativesYES)
 804     {
 805         nsets_dhdl = n_lambda_terms;
 806     }
 807     /* count the number of delta_g states */
 808     nsets_de = fep->lambda_stop_n - fep->lambda_start_n;
 809
 810     nsets = nsets_dhdl + nsets_de; /* dhdl + fep differences */
 811
 812     if (fep->n_lambda > 0 && (expand->elmcmove > elmcmoveNO))
 813     {
 814         nsets += 1;   /*add fep state for expanded ensemble */
 815     }
 816
 817     if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
 818     {
 819         nsets += 1;  /* add energy to the dhdl as well */
 820     }
 821
 822     nsetsextend = nsets;
 823     if ((ir->epc != epcNO) && (fep->n_lambda > 0) && (fep->init_lambda < 0))
 824     {
 825         nsetsextend += 1; /* for PV term, other terms possible if required for
 826                              the reduced potential (only needed with foreign
 827                              lambda, and only output when init_lambda is not
 828                              set in order to maintain compatibility of the
 829                              dhdl.xvg file) */
 830         write_pV     = true;
 831     }
 832     std::vector<std::string> setname(nsetsextend);
 833
 834     if (expand->elmcmove > elmcmoveNO)
 835     {
 836         /* state for the fep_vals, if we have alchemical sampling */
 837         setname[s++] = "Thermodynamic state";
 838     }
 839
 840     if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
 841     {
 842         std::string energy;
 843         switch (fep->edHdLPrintEnergy)
 844         {
 845             case edHdLPrintEnergyPOTENTIAL:
 846                 energy = gmx::formatString("%s (%s)", "Potential Energy", unit_energy);
 847                 break;
 848             case edHdLPrintEnergyTOTAL:
 849             case edHdLPrintEnergyYES:
 850             default:
 851                 energy = gmx::formatString("%s (%s)", "Total Energy", unit_energy);
 852         }
 853         setname[s++] = energy;
 854     }
 855
 856     if (fep->dhdl_derivatives == edhdlderivativesYES)
 857     {
 858         for (i = 0; i < efptNR; i++)
 859         {
 860             if (fep->separate_dvdl[i])
 861             {
 862                 std::string derivative;
 863                 if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))
 864                 {
 865                     /* compatibility output */
 866                     derivative = gmx::formatString("%s %s %.4f", dhdl, lambda, fep->init_lambda);
 867                 }
 868                 else
 869                 {
 870                     double lam = fep->init_lambda;
 871                     if (fep->init_lambda < 0)
 872                     {
 873                         lam = fep->all_lambda[i][fep->init_fep_state];
 874                     }
 875                     derivative = gmx::formatString("%s %s = %.4f", dhdl, efpt_singular_names[i],
 876                                                    lam);
 877                 }
 878                 setname[s++] = derivative;
 879             }
 880         }
 881     }
 882
 883     if (fep->n_lambda > 0)
 884     {
 885         /* g_bar has to determine the lambda values used in this simulation
 886          * from this xvg legend.
 887          */
 888
 889         if (expand->elmcmove > elmcmoveNO)
 890         {
 891             nsetsbegin = 1;  /* for including the expanded ensemble */
 892         }
 893         else
 894         {
 895             nsetsbegin = 0;
 896         }
 897
 898         if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
 899         {
 900             nsetsbegin += 1;
 901         }
 902         nsetsbegin += nsets_dhdl;
 903
 904         for (i = fep->lambda_start_n; i < fep->lambda_stop_n; i++)
 905         {
 906             print_lambda_vector(fep, i, false, false, lambda_vec_str);
 907             std::string buf;
 908             if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))
 909             {
 910                 /* for compatible dhdl.xvg files */
 911                 buf = gmx::formatString("%s %s %s", deltag, lambda, lambda_vec_str);
 912             }
 913             else
 914             {
 915                 buf = gmx::formatString("%s %s to %s", deltag, lambda, lambda_vec_str);
 916             }
 917
 918             if (ir->bSimTemp)
 919             {
 920                 /* print the temperature for this state if doing simulated annealing */
 921                 buf += gmx::formatString("T = %g (%s)",
 922                                          ir->simtempvals->temperatures[s-(nsetsbegin)],
 923                                          unit_temp_K);
 924             }
 925             setname[s++] = buf;
 926         }
 927         if (write_pV)
 928         {
 929             setname[s++] = gmx::formatString("pV (%s)", unit_energy);
 930         }
 931
 932         xvgrLegend(fp, setname, oenv);
 933     }
 934
 935     return fp;
 936 }
 937
 938 namespace gmx
 939 {
 940
 941 void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
 942                                        bool                    bSum,
 943                                        double                  time,
 944                                        real                    tmass,
 945                                        gmx_enerdata_t         *enerd,
 946                                        t_state                *state,
 947                                        t_lambda               *fep,
 948                                        t_expanded             *expand,
 949                                        matrix                  box,
 950                                        tensor                  svir,
 951                                        tensor                  fvir,
 952                                        tensor                  vir,
 953                                        tensor                  pres,
 954                                        gmx_ekindata_t         *ekind,
 955                                        rvec                    mu_tot,
 956                                        const gmx::Constraints *constr)
 957 {
 958     int    j, k, kk, n, gid;
 959     real   crmsd[2], tmp6[6];
 960     real   bs[tricl_boxs_nm.size()], vol, dens, pv, enthalpy;
 961     real   eee[egNR];
 962     double store_dhdl[efptNR];
 963     real   store_energy = 0;
 964     real   tmp;
 965
 966     /* Do NOT use the box in the state variable, but the separate box provided
 967      * as an argument. This is because we sometimes need to write the box from
 968      * the last timestep to match the trajectory frames.
 969      */
 970     add_ebin_indexed(ebin_, ie_, gmx::ArrayRef<bool>(bEner_), enerd->term, bSum);
 971     if (nCrmsd_)
 972     {
 973         crmsd[0] = constr->rmsd();
 974         add_ebin(ebin_, iconrmsd_, nCrmsd_, crmsd, false);
 975     }
 976     if (bDynBox_)
 977     {
 978         int nboxs;
 979         if (bTricl_)
 980         {
 981             bs[0] = box[XX][XX];
 982             bs[1] = box[YY][YY];
 983             bs[2] = box[ZZ][ZZ];
 984             bs[3] = box[YY][XX];
 985             bs[4] = box[ZZ][XX];
 986             bs[5] = box[ZZ][YY];
 987             nboxs = tricl_boxs_nm.size();
 988         }
 989         else
 990         {
 991             bs[0] = box[XX][XX];
 992             bs[1] = box[YY][YY];
 993             bs[2] = box[ZZ][ZZ];
 994             nboxs = boxs_nm.size();
 995         }
 996         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 997         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 998         add_ebin(ebin_, ib_, nboxs, bs, bSum);
 999         add_ebin(ebin_, ivol_, 1, &vol, bSum);
1000         add_ebin(ebin_, idens_, 1, &dens, bSum);
1001
1002         if (bDiagPres_)
1003         {
1004             /* This is pV (in kJ/mol).  The pressure is the reference pressure,
1005                not the instantaneous pressure */
1006             pv = vol*ref_p_/PRESFAC;
1007
1008             add_ebin(ebin_, ipv_, 1, &pv, bSum);
1009             enthalpy = pv + enerd->term[F_ETOT];
1010             add_ebin(ebin_, ienthalpy_, 1, &enthalpy, bSum);
1011         }
1012     }
1013     if (bConstrVir_)
1014     {
1015         add_ebin(ebin_, isvir_, 9, svir[0], bSum);
1016         add_ebin(ebin_, ifvir_, 9, fvir[0], bSum);
1017     }
1018     if (bPres_)
1019     {
1020         add_ebin(ebin_, ivir_, 9, vir[0], bSum);
1021         add_ebin(ebin_, ipres_, 9, pres[0], bSum);
1022         tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
1023         add_ebin(ebin_, isurft_, 1, &tmp, bSum);
1024     }
1025     if (epc_ == epcPARRINELLORAHMAN || epc_ == epcMTTK)
1026     {
1027         tmp6[0] = state->boxv[XX][XX];
1028         tmp6[1] = state->boxv[YY][YY];
1029         tmp6[2] = state->boxv[ZZ][ZZ];
1030         tmp6[3] = state->boxv[YY][XX];
1031         tmp6[4] = state->boxv[ZZ][XX];
1032         tmp6[5] = state->boxv[ZZ][YY];
1033         add_ebin(ebin_, ipc_, bTricl_ ? 6 : 3, tmp6, bSum);
1034     }
1035     if (bMu_)
1036     {
1037         add_ebin(ebin_, imu_, 3, mu_tot, bSum);
1038     }
1039     if (ekind && ekind->cosacc.cos_accel != 0)
1040     {
1041         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
1042         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
1043         add_ebin(ebin_, ivcos_, 1, &(ekind->cosacc.vcos), bSum);
1044         /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
1045         tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
1046                  *dens*gmx::square(box[ZZ][ZZ]*NANO/(2*M_PI)));
1047         add_ebin(ebin_, ivisc_, 1, &tmp, bSum);
1048     }
1049     if (nE_ > 1)
1050     {
1051         n = 0;
1052         for (int i = 0; (i < nEg_); i++)
1053         {
1054             for (j = i; (j < nEg_); j++)
1055             {
1056                 gid = GID(i, j, nEg_);
1057                 for (k = kk = 0; (k < egNR); k++)
1058                 {
1059                     if (bEInd_[k])
1060                     {
1061                         eee[kk++] = enerd->grpp.ener[k][gid];
1062                     }
1063                 }
1064                 add_ebin(ebin_, igrp_[n], nEc_, eee, bSum);
1065                 n++;
1066             }
1067         }
1068     }
1069
1070     if (ekind)
1071     {
1072         for (int i = 0; (i < nTC_); i++)
1073         {
1074             tmp_r_[i] = ekind->tcstat[i].T;
1075         }
1076         add_ebin(ebin_, itemp_, nTC_, tmp_r_, bSum);
1077
1078         if (etc_ == etcNOSEHOOVER)
1079         {
1080             /* whether to print Nose-Hoover chains: */
1081             if (bPrintNHChains_)
1082             {
1083                 if (bNHC_trotter_)
1084                 {
1085                     for (int i = 0; (i < nTC_); i++)
1086                     {
1087                         for (j = 0; j < nNHC_; j++)
1088                         {
1089                             k                  = i*nNHC_+j;
1090                             tmp_r_[2*k]        = state->nosehoover_xi[k];
1091                             tmp_r_[2*k+1]      = state->nosehoover_vxi[k];
1092                         }
1093                     }
1094                     add_ebin(ebin_, itc_, mde_n_, tmp_r_, bSum);
1095
1096                     if (bMTTK_)
1097                     {
1098                         for (int i = 0; (i < nTCP_); i++)
1099                         {
1100                             for (j = 0; j < nNHC_; j++)
1101                             {
1102                                 k                  = i*nNHC_+j;
1103                                 tmp_r_[2*k]        = state->nhpres_xi[k];
1104                                 tmp_r_[2*k+1]      = state->nhpres_vxi[k];
1105                             }
1106                         }
1107                         add_ebin(ebin_, itcb_, mdeb_n_, tmp_r_, bSum);
1108                     }
1109                 }
1110                 else
1111                 {
1112                     for (int i = 0; (i < nTC_); i++)
1113                     {
1114                         tmp_r_[2*i]   = state->nosehoover_xi[i];
1115                         tmp_r_[2*i+1] = state->nosehoover_vxi[i];
1116                     }
1117                     add_ebin(ebin_, itc_, mde_n_, tmp_r_, bSum);
1118                 }
1119             }
1120         }
1121         else if (etc_ == etcBERENDSEN || etc_ == etcYES ||
1122                  etc_ == etcVRESCALE)
1123         {
1124             for (int i = 0; (i < nTC_); i++)
1125             {
1126                 tmp_r_[i] = ekind->tcstat[i].lambda;
1127             }
1128             add_ebin(ebin_, itc_, nTC_, tmp_r_, bSum);
1129         }
1130     }
1131
1132     if (ekind && nU_ > 1)
1133     {
1134         for (int i = 0; (i < nU_); i++)
1135         {
1136             copy_rvec(ekind->grpstat[i].u, tmp_v_[i]);
1137         }
1138         add_ebin(ebin_, iu_, 3*nU_, tmp_v_[0], bSum);
1139     }
1140
1141     ebin_increase_count(1, ebin_, bSum);
1142
1143     // BAR + thermodynamic integration values
1144     if ((fp_dhdl_ || dhc_) && bDoDHDL)
1145     {
1146         for (gmx::index i = 0; i < static_cast<gmx::index>(enerd->enerpart_lambda.size()) - 1; i++)
1147         {
1148             /* zero for simulated tempering */
1149             dE_[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0];
1150             if (numTemperatures_ > 0)
1151             {
1152                 GMX_RELEASE_ASSERT(numTemperatures_ > state->fep_state, "Number of lambdas in state is bigger then in input record");
1153                 GMX_RELEASE_ASSERT(numTemperatures_ >= static_cast<gmx::index>(enerd->enerpart_lambda.size()) - 1, "Number of lambdas in energy data is bigger then in input record");
1154                 /* MRS: is this right, given the way we have defined the exchange probabilities? */
1155                 /* is this even useful to have at all? */
1156                 dE_[i] += (temperatures_[i]/
1157                            temperatures_[state->fep_state]-1.0)*
1158                     enerd->term[F_EKIN];
1159             }
1160         }
1161
1162         if (fp_dhdl_)
1163         {
1164             fprintf(fp_dhdl_, "%.4f", time);
1165             /* the current free energy state */
1166
1167             /* print the current state if we are doing expanded ensemble */
1168             if (expand->elmcmove > elmcmoveNO)
1169             {
1170                 fprintf(fp_dhdl_, " %4d", state->fep_state);
1171             }
1172             /* total energy (for if the temperature changes */
1173
1174             if (fep->edHdLPrintEnergy != edHdLPrintEnergyNO)
1175             {
1176                 switch (fep->edHdLPrintEnergy)
1177                 {
1178                     case edHdLPrintEnergyPOTENTIAL:
1179                         store_energy = enerd->term[F_EPOT];
1180                         break;
1181                     case edHdLPrintEnergyTOTAL:
1182                     case edHdLPrintEnergyYES:
1183                     default:
1184                         store_energy = enerd->term[F_ETOT];
1185                 }
1186                 fprintf(fp_dhdl_, " %#.8g", store_energy);
1187             }
1188
1189             if (fep->dhdl_derivatives == edhdlderivativesYES)
1190             {
1191                 for (int i = 0; i < efptNR; i++)
1192                 {
1193                     if (fep->separate_dvdl[i])
1194                     {
1195                         /* assumes F_DVDL is first */
1196                         fprintf(fp_dhdl_, " %#.8g", enerd->term[F_DVDL+i]);
1197                     }
1198                 }
1199             }
1200             for (int i = fep->lambda_start_n; i < fep->lambda_stop_n; i++)
1201             {
1202                 fprintf(fp_dhdl_, " %#.8g", dE_[i]);
1203             }
1204             if (bDynBox_ &&
1205                 bDiagPres_ &&
1206                 (epc_ != epcNO) &&
1207                 !enerd->enerpart_lambda.empty() &&
1208                 (fep->init_lambda < 0))
1209             {
1210                 fprintf(fp_dhdl_, " %#.8g", pv);  /* PV term only needed when
1211                                                           there are alternate state
1212                                                           lambda and we're not in
1213                                                           compatibility mode */
1214             }
1215             fprintf(fp_dhdl_, "\n");
1216             /* and the binary free energy output */
1217         }
1218         if (dhc_ && bDoDHDL)
1219         {
1220             int idhdl = 0;
1221             for (int i = 0; i < efptNR; i++)
1222             {
1223                 if (fep->separate_dvdl[i])
1224                 {
1225                     /* assumes F_DVDL is first */
1226                     store_dhdl[idhdl] = enerd->term[F_DVDL+i];
1227                     idhdl            += 1;
1228                 }
1229             }
1230             store_energy = enerd->term[F_ETOT];
1231             /* store_dh is dE */
1232             mde_delta_h_coll_add_dh(dhc_,
1233                                     static_cast<double>(state->fep_state),
1234                                     store_energy,
1235                                     pv,
1236                                     store_dhdl,
1237                                     dE_ + fep->lambda_start_n,
1238                                     time);
1239         }
1240     }
1241     ;
1242 }
1243
1244 void EnergyOutput::recordNonEnergyStep()
1245 {
1246     ebin_increase_count(1, ebin_, false);
1247 }
1248
1249 void EnergyOutput::printHeader(FILE *log, int64_t steps, double time)
1250 {
1251     char buf[22];
1252
1253     fprintf(log, "   %12s   %12s\n"
1254             "   %12s   %12.5f\n\n",
1255             "Step", "Time", gmx_step_str(steps, buf), time);
1256 }
1257
1258 void EnergyOutput::printStepToEnergyFile(ener_file *fp_ene, bool bEne, bool bDR, bool bOR,
1259                                          FILE *log,
1260                                          int64_t step, double time,
1261                                          t_fcdata *fcd,
1262                                          gmx::Awh *awh)
1263 {
1264     t_enxframe   fr;
1265     init_enxframe(&fr);
1266     fr.t            = time;
1267     fr.step         = step;
1268     fr.nsteps       = ebin_->nsteps;
1269     fr.dt           = delta_t_;
1270     fr.nsum         = ebin_->nsum;
1271     fr.nre          = (bEne) ? ebin_->nener : 0;
1272     fr.ener         = ebin_->e;
1273     int ndisre          = bDR ? fcd->disres.npair : 0;
1274     /* these are for the old-style blocks (1 subblock, only reals), because
1275        there can be only one per ID for these */
1276     int          nr[enxNR];
1277     int          id[enxNR];
1278     real        *block[enxNR];
1279     /* Optional additional old-style (real-only) blocks. */
1280     for (int i = 0; i < enxNR; i++)
1281     {
1282         nr[i] = 0;
1283     }
1284
1285     if (bOR && fcd->orires.nr > 0)
1286     {
1287         diagonalize_orires_tensors(&(fcd->orires));
1288         nr[enxOR]     = fcd->orires.nr;
1289         block[enxOR]  = fcd->orires.otav;
1290         id[enxOR]     = enxOR;
1291         nr[enxORI]    = (fcd->orires.oinsl != fcd->orires.otav) ?
1292             fcd->orires.nr : 0;
1293         block[enxORI] = fcd->orires.oinsl;
1294         id[enxORI]    = enxORI;
1295         nr[enxORT]    = fcd->orires.nex*12;
1296         block[enxORT] = fcd->orires.eig;
1297         id[enxORT]    = enxORT;
1298     }
1299
1300     /* whether we are going to write anything out: */
1301     if (fr.nre || ndisre || nr[enxOR] || nr[enxORI])
1302     {
1303         /* the old-style blocks go first */
1304         fr.nblock = 0;
1305         for (int i = 0; i < enxNR; i++)
1306         {
1307             if (nr[i] > 0)
1308             {
1309                 fr.nblock = i + 1;
1310             }
1311         }
1312         add_blocks_enxframe(&fr, fr.nblock);
1313         for (int b = 0; b < fr.nblock; b++)
1314         {
1315             add_subblocks_enxblock(&(fr.block[b]), 1);
1316             fr.block[b].id        = id[b];
1317             fr.block[b].sub[0].nr = nr[b];
1318 #if !GMX_DOUBLE
1319             fr.block[b].sub[0].type = xdr_datatype_float;
1320             fr.block[b].sub[0].fval = block[b];
1321 #else
1322             fr.block[b].sub[0].type = xdr_datatype_double;
1323             fr.block[b].sub[0].dval = block[b];
1324 #endif
1325         }
1326
1327         /* check for disre block & fill it. */
1328         if (ndisre > 0)
1329         {
1330             int db = fr.nblock;
1331             fr.nblock += 1;
1332             add_blocks_enxframe(&fr, fr.nblock);
1333
1334             add_subblocks_enxblock(&(fr.block[db]), 2);
1335             fr.block[db].id        = enxDISRE;
1336             fr.block[db].sub[0].nr = ndisre;
1337             fr.block[db].sub[1].nr = ndisre;
1338 #if !GMX_DOUBLE
1339             fr.block[db].sub[0].type = xdr_datatype_float;
1340             fr.block[db].sub[1].type = xdr_datatype_float;
1341             fr.block[db].sub[0].fval = fcd->disres.rt;
1342             fr.block[db].sub[1].fval = fcd->disres.rm3tav;
1343 #else
1344             fr.block[db].sub[0].type = xdr_datatype_double;
1345             fr.block[db].sub[1].type = xdr_datatype_double;
1346             fr.block[db].sub[0].dval = fcd->disres.rt;
1347             fr.block[db].sub[1].dval = fcd->disres.rm3tav;
1348 #endif
1349         }
1350         /* here we can put new-style blocks */
1351
1352         /* Free energy perturbation blocks */
1353         if (dhc_)
1354         {
1355             mde_delta_h_coll_handle_block(dhc_, &fr, fr.nblock);
1356         }
1357
1358         /* we can now free & reset the data in the blocks */
1359         if (dhc_)
1360         {
1361             mde_delta_h_coll_reset(dhc_);
1362         }
1363
1364         /* AWH bias blocks. */
1365         if (awh != nullptr)  // TODO: add boolean flag.
1366         {
1367             awh->writeToEnergyFrame(step, &fr);
1368         }
1369
1370         /* do the actual I/O */
1371         do_enx(fp_ene, &fr);
1372         if (fr.nre)
1373         {
1374             /* We have stored the sums, so reset the sum history */
1375             reset_ebin_sums(ebin_);
1376         }
1377     }
1378     free_enxframe(&fr);
1379     if (log)
1380     {
1381         if (bOR && fcd->orires.nr > 0)
1382         {
1383             print_orires_log(log, &(fcd->orires));
1384         }
1385
1386         fprintf(log, "   Energies (%s)\n", unit_energy);
1387         pr_ebin(log, ebin_, ie_, f_nre_+nCrmsd_, 5, eprNORMAL, true);
1388         fprintf(log, "\n");
1389     }
1390 }
1391
1392 void EnergyOutput::printAnnealingTemperatures(FILE *log, SimulationGroups *groups, t_grpopts *opts)
1393 {
1394     if (log)
1395     {
1396         if (opts)
1397         {
1398             for (int i = 0; i < opts->ngtc; i++)
1399             {
1400                 if (opts->annealing[i] != eannNO)
1401                 {
1402                     fprintf(log, "Current ref_t for group %s: %8.1f\n",
1403                             *(groups->groupNames[groups->groups[SimulationAtomGroupType::TemperatureCoupling][i]]),
1404                             opts->ref_t[i]);
1405                 }
1406             }
1407             fprintf(log, "\n");
1408         }
1409     }
1410 }
1411
1412 void EnergyOutput::printAverages(FILE *log, SimulationGroups *groups)
1413 {
1414     if (ebin_->nsum_sim <= 0)
1415     {
1416         if (log)
1417         {
1418             fprintf(log, "Not enough data recorded to report energy averages\n");
1419         }
1420         return;
1421     }
1422     if (log)
1423     {
1424
1425         char buf1[22], buf2[22];
1426
1427         fprintf(log, "\t<======  ###############  ==>\n");
1428         fprintf(log, "\t<====  A V E R A G E S  ====>\n");
1429         fprintf(log, "\t<==  ###############  ======>\n\n");
1430
1431         fprintf(log, "\tStatistics over %s steps using %s frames\n",
1432                 gmx_step_str(ebin_->nsteps_sim, buf1),
1433                 gmx_step_str(ebin_->nsum_sim, buf2));
1434         fprintf(log, "\n");
1435
1436         fprintf(log, "   Energies (%s)\n", unit_energy);
1437         pr_ebin(log, ebin_, ie_, f_nre_+nCrmsd_, 5, eprAVER, true);
1438         fprintf(log, "\n");
1439
1440         if (bDynBox_)
1441         {
1442             pr_ebin(log, ebin_, ib_, bTricl_ ? tricl_boxs_nm.size() : boxs_nm.size(), 5,
1443                     eprAVER, true);
1444             fprintf(log, "\n");
1445         }
1446         if (bConstrVir_)
1447         {
1448             fprintf(log, "   Constraint Virial (%s)\n", unit_energy);
1449             pr_ebin(log, ebin_, isvir_, 9, 3, eprAVER, false);
1450             fprintf(log, "\n");
1451             fprintf(log, "   Force Virial (%s)\n", unit_energy);
1452             pr_ebin(log, ebin_, ifvir_, 9, 3, eprAVER, false);
1453             fprintf(log, "\n");
1454         }
1455         if (bPres_)
1456         {
1457             fprintf(log, "   Total Virial (%s)\n", unit_energy);
1458             pr_ebin(log, ebin_, ivir_, 9, 3, eprAVER, false);
1459             fprintf(log, "\n");
1460             fprintf(log, "   Pressure (%s)\n", unit_pres_bar);
1461             pr_ebin(log, ebin_, ipres_, 9, 3, eprAVER, false);
1462             fprintf(log, "\n");
1463         }
1464         if (bMu_)
1465         {
1466             fprintf(log, "   Total Dipole (%s)\n", unit_dipole_D);
1467             pr_ebin(log, ebin_, imu_, 3, 3, eprAVER, false);
1468             fprintf(log, "\n");
1469         }
1470
1471         if (nE_ > 1)
1472         {
1473             int padding = 8 - strlen(unit_energy);
1474             fprintf(log, "%*sEpot (%s)   ", padding, "", unit_energy);
1475             for (int i = 0; (i < egNR); i++)
1476             {
1477                 if (bEInd_[i])
1478                 {
1479                     fprintf(log, "%12s   ", egrp_nm[i]);
1480                 }
1481             }
1482             fprintf(log, "\n");
1483
1484             int n = 0;
1485             for (int i = 0; (i < nEg_); i++)
1486             {
1487                 int ni = groups->groups[SimulationAtomGroupType::EnergyOutput][i];
1488                 for (int j = i; (j < nEg_); j++)
1489                 {
1490                     int nj      = groups->groups[SimulationAtomGroupType::EnergyOutput][j];
1491                     int padding = 14 - (strlen(*(groups->groupNames[ni])) +
1492                                         strlen(*(groups->groupNames[nj])));
1493                     fprintf(log, "%*s%s-%s", padding, "",
1494                             *(groups->groupNames[ni]),
1495                             *(groups->groupNames[nj]));
1496                     pr_ebin(log, ebin_, igrp_[n], nEc_, nEc_, eprAVER,
1497                             false);
1498                     n++;
1499                 }
1500             }
1501             fprintf(log, "\n");
1502         }
1503         if (nTC_ > 1)
1504         {
1505             pr_ebin(log, ebin_, itemp_, nTC_, 4, eprAVER, true);
1506             fprintf(log, "\n");
1507         }
1508         if (nU_ > 1)
1509         {
1510             fprintf(log, "%15s   %12s   %12s   %12s\n",
1511                     "Group", "Ux", "Uy", "Uz");
1512             for (int i = 0; (i < nU_); i++)
1513             {
1514                 int ni = groups->groups[SimulationAtomGroupType::Acceleration][i];
1515                 fprintf(log, "%15s", *groups->groupNames[ni]);
1516                 pr_ebin(log, ebin_, iu_+3*i, 3, 3, eprAVER, false);
1517             }
1518             fprintf(log, "\n");
1519         }
1520     }
1521 }
1522
1523 void EnergyOutput::fillEnergyHistory(energyhistory_t *enerhist) const
1524 {
1525     const t_ebin * const ebin = ebin_;
1526
1527     enerhist->nsteps     = ebin->nsteps;
1528     enerhist->nsum       = ebin->nsum;
1529     enerhist->nsteps_sim = ebin->nsteps_sim;
1530     enerhist->nsum_sim   = ebin->nsum_sim;
1531
1532     if (ebin->nsum > 0)
1533     {
1534         /* This will only actually resize the first time */
1535         enerhist->ener_ave.resize(ebin->nener);
1536         enerhist->ener_sum.resize(ebin->nener);
1537
1538         for (int i = 0; i < ebin->nener; i++)
1539         {
1540             enerhist->ener_ave[i] = ebin->e[i].eav;
1541             enerhist->ener_sum[i] = ebin->e[i].esum;
1542         }
1543     }
1544
1545     if (ebin->nsum_sim > 0)
1546     {
1547         /* This will only actually resize the first time */
1548         enerhist->ener_sum_sim.resize(ebin->nener);
1549
1550         for (int i = 0; i < ebin->nener; i++)
1551         {
1552             enerhist->ener_sum_sim[i] = ebin->e_sim[i].esum;
1553         }
1554     }
1555     if (dhc_)
1556     {
1557         mde_delta_h_coll_update_energyhistory(dhc_, enerhist);
1558     }
1559 }
1560
1561 void EnergyOutput::restoreFromEnergyHistory(const energyhistory_t &enerhist)
1562 {
1563     unsigned int nener = static_cast<unsigned int>(ebin_->nener);
1564
1565     if ((enerhist.nsum     > 0 && nener != enerhist.ener_sum.size()) ||
1566         (enerhist.nsum_sim > 0 && nener != enerhist.ener_sum_sim.size()))
1567     {
1568         gmx_fatal(FARGS, "Mismatch between number of energies in run input (%u) and checkpoint file (%zu or %zu).",
1569                   nener, enerhist.ener_sum.size(), enerhist.ener_sum_sim.size());
1570     }
1571
1572     ebin_->nsteps     = enerhist.nsteps;
1573     ebin_->nsum       = enerhist.nsum;
1574     ebin_->nsteps_sim = enerhist.nsteps_sim;
1575     ebin_->nsum_sim   = enerhist.nsum_sim;
1576
1577     for (int i = 0; i < ebin_->nener; i++)
1578     {
1579         ebin_->e[i].eav  =
1580             (enerhist.nsum > 0 ? enerhist.ener_ave[i] : 0);
1581         ebin_->e[i].esum =
1582             (enerhist.nsum > 0 ? enerhist.ener_sum[i] : 0);
1583         ebin_->e_sim[i].esum =
1584             (enerhist.nsum_sim > 0 ? enerhist.ener_sum_sim[i] : 0);
1585     }
1586     if (dhc_)
1587     {
1588         mde_delta_h_coll_restore_energyhistory(dhc_, enerhist.deltaHForeignLambdas.get());
1589     }
1590 }
1591
1592 int EnergyOutput::numEnergyTerms() const
1593 {
1594     return ebin_->nener;
1595 }
1596
1597 } // namespace gmx