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