src/gromacs/mdlib/mdebin.c

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