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