src/gromacs/mdlib/mdebin.c

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