src/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 "mdrun.h"
  59 #include "mdebin_bar.h"
  60
  61
  62 static const char *conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
  63
  64 static const char *boxs_nm[] = { "Box-X", "Box-Y", "Box-Z" };
  65
  66 static const char *tricl_boxs_nm[] = {
  67     "Box-XX", "Box-YY", "Box-ZZ",
  68     "Box-YX", "Box-ZX", "Box-ZY"
  69 };
  70
  71 static const char *vol_nm[] = { "Volume" };
  72
  73 static const char *dens_nm[] = {"Density" };
  74
  75 static const char *pv_nm[] = {"pV" };
  76
  77 static const char *enthalpy_nm[] = {"Enthalpy" };
  78
  79 static const char *boxvel_nm[] = {
  80     "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
  81     "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
  82 };
  83
  84 #define NBOXS asize(boxs_nm)
  85 #define NTRICLBOXS asize(tricl_boxs_nm)
  86
  87 t_mdebin *init_mdebin(ener_file_t fp_ene,
  88                       const gmx_mtop_t *mtop,
  89                       const t_inputrec *ir,
  90                       FILE *fp_dhdl)
  91 {
  92     const char *ener_nm[F_NRE];
  93     static const char *vir_nm[] = {
  94         "Vir-XX", "Vir-XY", "Vir-XZ",
  95         "Vir-YX", "Vir-YY", "Vir-YZ",
  96         "Vir-ZX", "Vir-ZY", "Vir-ZZ"
  97     };
  98     static const char *sv_nm[] = {
  99         "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
 100         "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
 101         "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
 102     };
 103     static const char *fv_nm[] = {
 104         "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
 105         "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
 106         "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
 107     };
 108     static const char *pres_nm[] = {
 109         "Pres-XX","Pres-XY","Pres-XZ",
 110         "Pres-YX","Pres-YY","Pres-YZ",
 111         "Pres-ZX","Pres-ZY","Pres-ZZ"
 112     };
 113     static const char *surft_nm[] = {
 114         "#Surf*SurfTen"
 115     };
 116     static const char *mu_nm[] = {
 117         "Mu-X", "Mu-Y", "Mu-Z"
 118     };
 119     static const char *vcos_nm[] = {
 120         "2CosZ*Vel-X"
 121     };
 122     static const char *visc_nm[] = {
 123         "1/Viscosity"
 124     };
 125     static const char *baro_nm[] = {
 126         "Barostat"
 127     };
 128
 129     char     **grpnms;
 130     const gmx_groups_t *groups;
 131     char     **gnm;
 132     char     buf[256];
 133     const char     *bufi;
 134     t_mdebin *md;
 135     int      i,j,ni,nj,n,nh,k,kk,ncon,nset;
 136     gmx_bool     bBHAM,bNoseHoover,b14;
 137
 138     snew(md,1);
 139
 140     md->bVir=TRUE;
 141     md->bPress=TRUE;
 142     md->bSurft=TRUE;
 143     md->bMu=TRUE;
 144
 145     if (EI_DYNAMICS(ir->eI))
 146     {
 147         md->delta_t = ir->delta_t;
 148     }
 149     else
 150     {
 151         md->delta_t = 0;
 152     }
 153
 154     groups = &mtop->groups;
 155
 156     bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
 157     b14   = (gmx_mtop_ftype_count(mtop,F_LJ14) > 0 ||
 158              gmx_mtop_ftype_count(mtop,F_LJC14_Q) > 0);
 159
 160     ncon = gmx_mtop_ftype_count(mtop,F_CONSTR);
 161     nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
 162     md->bConstr    = (ncon > 0 || nset > 0);
 163     md->bConstrVir = FALSE;
 164     if (md->bConstr) {
 165         if (ncon > 0 && ir->eConstrAlg == econtLINCS) {
 166             if (ir->eI == eiSD2)
 167                 md->nCrmsd = 2;
 168             else
 169                 md->nCrmsd = 1;
 170         }
 171         md->bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL);
 172     } else {
 173         md->nCrmsd = 0;
 174     }
 175
 176     /* Energy monitoring */
 177     for(i=0;i<egNR;i++)
 178     {
 179         md->bEInd[i]=FALSE;
 180     }
 181
 182 #ifndef GMX_OPENMM
 183     for(i=0; i<F_NRE; i++)
 184     {
 185         md->bEner[i] = FALSE;
 186         if (i == F_LJ)
 187             md->bEner[i] = !bBHAM;
 188         else if (i == F_BHAM)
 189             md->bEner[i] = bBHAM;
 190         else if (i == F_EQM)
 191             md->bEner[i] = ir->bQMMM;
 192         else if (i == F_COUL_LR)
 193             md->bEner[i] = (ir->rcoulomb > ir->rlist);
 194         else if (i == F_LJ_LR)
 195             md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist);
 196         else if (i == F_BHAM_LR)
 197             md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist);
 198         else if (i == F_RF_EXCL)
 199             md->bEner[i] = (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC);
 200         else if (i == F_COUL_RECIP)
 201             md->bEner[i] = EEL_FULL(ir->coulombtype);
 202         else if (i == F_LJ14)
 203             md->bEner[i] = b14;
 204         else if (i == F_COUL14)
 205             md->bEner[i] = b14;
 206         else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
 207             md->bEner[i] = FALSE;
 208         else if ((i == F_DVDL_COUL && ir->fepvals->separate_dvdl[efptCOUL]) ||
 209                  (i == F_DVDL_VDW  && ir->fepvals->separate_dvdl[efptVDW]) ||
 210                  (i == F_DVDL_BONDED && ir->fepvals->separate_dvdl[efptBONDED]) ||
 211                  (i == F_DVDL_RESTRAINT && ir->fepvals->separate_dvdl[efptRESTRAINT]) ||
 212                  (i == F_DKDL && ir->fepvals->separate_dvdl[efptMASS]) ||
 213                  (i == F_DVDL && ir->fepvals->separate_dvdl[efptFEP]))
 214             md->bEner[i] = (ir->efep != efepNO);
 215         else if ((interaction_function[i].flags & IF_VSITE) ||
 216                  (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
 217             md->bEner[i] = FALSE;
 218         else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES)  || (i==F_EQM))
 219             md->bEner[i] = TRUE;
 220         else if ((i == F_GBPOL) && ir->implicit_solvent==eisGBSA)
 221             md->bEner[i] = TRUE;
 222         else if ((i == F_NPSOLVATION) && ir->implicit_solvent==eisGBSA && (ir->sa_algorithm != esaNO))
 223             md->bEner[i] = TRUE;
 224         else if ((i == F_GB12) || (i == F_GB13) || (i == F_GB14))
 225             md->bEner[i] = FALSE;
 226         else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
 227             md->bEner[i] = EI_DYNAMICS(ir->eI);
 228         else if (i==F_VTEMP)
 229             md->bEner[i] =  (EI_DYNAMICS(ir->eI) && getenv("GMX_VIRIAL_TEMPERATURE"));
 230         else if (i == F_DISPCORR || i == F_PDISPCORR)
 231             md->bEner[i] = (ir->eDispCorr != edispcNO);
 232         else if (i == F_DISRESVIOL)
 233             md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_DISRES) > 0);
 234         else if (i == F_ORIRESDEV)
 235             md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0);
 236         else if (i == F_CONNBONDS)
 237             md->bEner[i] = FALSE;
 238         else if (i == F_COM_PULL)
 239             md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F || ir->bRot);
 240         else if (i == F_ECONSERVED)
 241             md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) &&
 242                             (ir->epc == epcNO || ir->epc==epcMTTK));
 243         else
 244             md->bEner[i] = (gmx_mtop_ftype_count(mtop,i) > 0);
 245     }
 246 #else
 247     /* OpenMM always produces only the following 4 energy terms */
 248     md->bEner[F_EPOT] = TRUE;
 249     md->bEner[F_EKIN] = TRUE;
 250     md->bEner[F_ETOT] = TRUE;
 251     md->bEner[F_TEMP] = TRUE;
 252 #endif
 253
 254     /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
 255     if (ir->bAdress && !debug) {
 256         for (i = 0; i < F_NRE; i++) {
 257             md->bEner[i] = FALSE;
 258             if(i == F_EKIN){ md->bEner[i] = TRUE;}
 259             if(i == F_TEMP){ md->bEner[i] = TRUE;}
 260         }
 261         md->bVir=FALSE;
 262         md->bPress=FALSE;
 263         md->bSurft=FALSE;
 264         md->bMu=FALSE;
 265     }
 266
 267     md->f_nre=0;
 268     for(i=0; i<F_NRE; i++)
 269     {
 270         if (md->bEner[i])
 271         {
 272             /* FIXME: The constness should not be cast away */
 273             /*ener_nm[f_nre]=(char *)interaction_function[i].longname;*/
 274             ener_nm[md->f_nre]=interaction_function[i].longname;
 275             md->f_nre++;
 276         }
 277     }
 278
 279     md->epc = ir->epc;
 280     md->bDiagPres = !TRICLINIC(ir->ref_p);
 281     md->ref_p = (ir->ref_p[XX][XX]+ir->ref_p[YY][YY]+ir->ref_p[ZZ][ZZ])/DIM;
 282     md->bTricl = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
 283     md->bDynBox = DYNAMIC_BOX(*ir);
 284     md->etc = ir->etc;
 285     md->bNHC_trotter = IR_NVT_TROTTER(ir);
 286     md->bPrintNHChains = ir-> bPrintNHChains;
 287     md->bMTTK = (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(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         md->imu    = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
 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         add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
 885     if (ekind && ekind->cosacc.cos_accel != 0)
 886     {
 887         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 888         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 889         add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
 890         /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
 891         tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
 892                  *dens*vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
 893         add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
 894     }
 895     if (md->nE > 1)
 896     {
 897         n=0;
 898         for(i=0; (i<md->nEg); i++)
 899         {
 900             for(j=i; (j<md->nEg); j++)
 901             {
 902                 gid=GID(i,j,md->nEg);
 903                 for(k=kk=0; (k<egNR); k++)
 904                 {
 905                     if (md->bEInd[k])
 906                     {
 907                         eee[kk++] = enerd->grpp.ener[k][gid];
 908                     }
 909                 }
 910                 add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
 911                 n++;
 912             }
 913         }
 914     }
 915
 916     if (ekind)
 917     {
 918         for(i=0; (i<md->nTC); i++)
 919         {
 920             md->tmp_r[i] = ekind->tcstat[i].T;
 921         }
 922         add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
 923
 924         if (md->etc == etcNOSEHOOVER)
 925         {
 926             /* whether to print Nose-Hoover chains: */
 927             if (md->bPrintNHChains)
 928             {
 929                 if (md->bNHC_trotter)
 930                 {
 931                     for(i=0; (i<md->nTC); i++)
 932                     {
 933                         for (j=0;j<md->nNHC;j++)
 934                         {
 935                             k = i*md->nNHC+j;
 936                             md->tmp_r[2*k] = state->nosehoover_xi[k];
 937                             md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
 938                         }
 939                     }
 940                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
 941
 942                     if (md->bMTTK) {
 943                         for(i=0; (i<md->nTCP); i++)
 944                         {
 945                             for (j=0;j<md->nNHC;j++)
 946                             {
 947                                 k = i*md->nNHC+j;
 948                                 md->tmp_r[2*k] = state->nhpres_xi[k];
 949                                 md->tmp_r[2*k+1] = state->nhpres_vxi[k];
 950                             }
 951                         }
 952                         add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
 953                     }
 954                 }
 955                 else
 956                 {
 957                     for(i=0; (i<md->nTC); i++)
 958                     {
 959                         md->tmp_r[2*i] = state->nosehoover_xi[i];
 960                         md->tmp_r[2*i+1] = state->nosehoover_vxi[i];
 961                     }
 962                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
 963                 }
 964             }
 965         }
 966         else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
 967                  md->etc == etcVRESCALE)
 968         {
 969             for(i=0; (i<md->nTC); i++)
 970             {
 971                 md->tmp_r[i] = ekind->tcstat[i].lambda;
 972             }
 973             add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
 974         }
 975     }
 976
 977     if (ekind && md->nU > 1)
 978     {
 979         for(i=0; (i<md->nU); i++)
 980         {
 981             copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
 982         }
 983         add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
 984     }
 985
 986     ebin_increase_count(md->ebin,bSum);
 987
 988     /* BAR + thermodynamic integration values */
 989     if ((md->fp_dhdl || md->dhc) && bDoDHDL && (enerd->n_lambda > 0))
 990     {
 991         snew(dE,enerd->n_lambda-1);
 992         for(i=0; i<enerd->n_lambda-1; i++) {
 993             dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0];  /* zero for simulated tempering */
 994             if (md->temperatures!=NULL)
 995             {
 996                 /* MRS: is this right, given the way we have defined the exchange probabilities? */
 997                 /* is this even useful to have at all? */
 998                 dE[i] += (md->temperatures[i]/md->temperatures[state->fep_state]-1.0)*enerd->term[F_EKIN];
 999             }
1000         }
1001     }
1002
1003     if (md->fp_dhdl && bDoDHDL)
1004     {
1005         fprintf(md->fp_dhdl,"%.4f",time);
1006         /* the current free energy state */
1007
1008         /* print the current state if we are doing expanded ensemble */
1009         if (expand->elmcmove > elmcmoveNO) {
1010             fprintf(md->fp_dhdl," %4d",state->fep_state);
1011         }
1012         /* total energy (for if the temperature changes */
1013         if (fep->bPrintEnergy)
1014         {
1015             store_energy = enerd->term[F_ETOT];
1016             fprintf(md->fp_dhdl," %#.8g",store_energy);
1017         }
1018
1019         for (i=0;i<efptNR;i++)
1020         {
1021             if (fep->separate_dvdl[i])
1022             {
1023                 fprintf(md->fp_dhdl," %#.8g",enerd->term[F_DVDL+i]); /* assumes F_DVDL is first */
1024             }
1025         }
1026         for(i=1; i<enerd->n_lambda; i++)
1027         {
1028             fprintf(md->fp_dhdl," %#.8g",dE[i-1]);
1029
1030         }
1031         if ((md->epc!=epcNO)  && (enerd->n_lambda > 0))
1032         {
1033             fprintf(md->fp_dhdl," %#.8g",pv);   /* PV term only needed when there are alternate state lambda */
1034         }
1035         fprintf(md->fp_dhdl,"\n");
1036         /* and the binary free energy output */
1037     }
1038     if (md->dhc && bDoDHDL)
1039     {
1040         int idhdl = 0;
1041         for (i=0;i<efptNR;i++)
1042         {
1043             if (fep->separate_dvdl[i])
1044             {
1045                 store_dhdl[idhdl] = enerd->term[F_DVDL+i]; /* assumes F_DVDL is first */
1046                 idhdl+=1;
1047             }
1048         }
1049         /* store_dh is dE */
1050         mde_delta_h_coll_add_dh(md->dhc,
1051                                 (double)state->fep_state,
1052                                 store_energy,
1053                                 pv,
1054                                 (expand->elamstats>elamstatsNO),
1055                                 (fep->bPrintEnergy),
1056                                 (md->epc!=epcNO),
1057                                 idhdl,
1058                                 fep->n_lambda,
1059                                 store_dhdl,
1060                                 dE,
1061                                 time);
1062     }
1063     if ((md->fp_dhdl || md->dhc) && bDoDHDL && (enerd->n_lambda >0))
1064     {
1065         sfree(dE);
1066     }
1067 }
1068
1069
1070 void upd_mdebin_step(t_mdebin *md)
1071 {
1072     ebin_increase_count(md->ebin,FALSE);
1073 }
1074
1075 static void npr(FILE *log,int n,char c)
1076 {
1077     for(; (n>0); n--) fprintf(log,"%c",c);
1078 }
1079
1080 static void pprint(FILE *log,const char *s,t_mdebin *md)
1081 {
1082     char CHAR='#';
1083     int  slen;
1084     char buf1[22],buf2[22];
1085
1086     slen = strlen(s);
1087     fprintf(log,"\t<======  ");
1088     npr(log,slen,CHAR);
1089     fprintf(log,"  ==>\n");
1090     fprintf(log,"\t<====  %s  ====>\n",s);
1091     fprintf(log,"\t<==  ");
1092     npr(log,slen,CHAR);
1093     fprintf(log,"  ======>\n\n");
1094
1095     fprintf(log,"\tStatistics over %s steps using %s frames\n",
1096             gmx_step_str(md->ebin->nsteps_sim,buf1),
1097             gmx_step_str(md->ebin->nsum_sim,buf2));
1098     fprintf(log,"\n");
1099 }
1100
1101 void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lambda)
1102 {
1103     char buf[22];
1104
1105     fprintf(log,"   %12s   %12s   %12s\n"
1106             "   %12s   %12.5f   %12.5f\n\n",
1107             "Step","Time","Lambda",gmx_step_str(steps,buf),time,lambda);
1108 }
1109
1110 void print_ebin(ener_file_t fp_ene,gmx_bool bEne,gmx_bool bDR,gmx_bool bOR,
1111                 FILE *log,
1112                 gmx_large_int_t step,double time,
1113                 int mode,gmx_bool bCompact,
1114                 t_mdebin *md,t_fcdata *fcd,
1115                 gmx_groups_t *groups,t_grpopts *opts)
1116 {
1117     /*static char **grpnms=NULL;*/
1118     char        buf[246];
1119     int         i,j,n,ni,nj,ndr,nor,b;
1120     int         ndisre=0;
1121     real        *disre_rm3tav, *disre_rt;
1122
1123     /* these are for the old-style blocks (1 subblock, only reals), because
1124        there can be only one per ID for these */
1125     int         nr[enxNR];
1126     int         id[enxNR];
1127     real        *block[enxNR];
1128
1129     /* temporary arrays for the lambda values to write out */
1130     double      enxlambda_data[2];
1131
1132     t_enxframe  fr;
1133
1134     switch (mode)
1135     {
1136         case eprNORMAL:
1137             init_enxframe(&fr);
1138             fr.t            = time;
1139             fr.step         = step;
1140             fr.nsteps       = md->ebin->nsteps;
1141             fr.dt           = md->delta_t;
1142             fr.nsum         = md->ebin->nsum;
1143             fr.nre          = (bEne) ? md->ebin->nener : 0;
1144             fr.ener         = md->ebin->e;
1145             ndisre          = bDR ? fcd->disres.npair : 0;
1146             disre_rm3tav    = fcd->disres.rm3tav;
1147             disre_rt        = fcd->disres.rt;
1148             /* Optional additional old-style (real-only) blocks. */
1149             for(i=0; i<enxNR; i++)
1150             {
1151                 nr[i] = 0;
1152             }
1153             if (fcd->orires.nr > 0 && bOR)
1154             {
1155                 diagonalize_orires_tensors(&(fcd->orires));
1156                 nr[enxOR]     = fcd->orires.nr;
1157                 block[enxOR]  = fcd->orires.otav;
1158                 id[enxOR]     = enxOR;
1159                 nr[enxORI]    = (fcd->orires.oinsl != fcd->orires.otav) ?
1160                           fcd->orires.nr : 0;
1161                 block[enxORI] = fcd->orires.oinsl;
1162                 id[enxORI]    = enxORI;
1163                 nr[enxORT]    = fcd->orires.nex*12;
1164                 block[enxORT] = fcd->orires.eig;
1165                 id[enxORT]    = enxORT;
1166             }
1167
1168             /* whether we are going to wrte anything out: */
1169             if (fr.nre || ndisre || nr[enxOR] || nr[enxORI])
1170             {
1171
1172                 /* the old-style blocks go first */
1173                 fr.nblock = 0;
1174                 for(i=0; i<enxNR; i++)
1175                 {
1176                     if (nr[i] > 0)
1177                     {
1178                         fr.nblock = i + 1;
1179                     }
1180                 }
1181                 add_blocks_enxframe(&fr, fr.nblock);
1182                 for(b=0;b<fr.nblock;b++)
1183                 {
1184                     add_subblocks_enxblock(&(fr.block[b]), 1);
1185                     fr.block[b].id=id[b];
1186                     fr.block[b].sub[0].nr = nr[b];
1187 #ifndef GMX_DOUBLE
1188                     fr.block[b].sub[0].type = xdr_datatype_float;
1189                     fr.block[b].sub[0].fval = block[b];
1190 #else
1191                     fr.block[b].sub[0].type = xdr_datatype_double;
1192                     fr.block[b].sub[0].dval = block[b];
1193 #endif
1194                 }
1195
1196                 /* check for disre block & fill it. */
1197                 if (ndisre>0)
1198                 {
1199                     int db = fr.nblock;
1200                     fr.nblock+=1;
1201                     add_blocks_enxframe(&fr, fr.nblock);
1202
1203                     add_subblocks_enxblock(&(fr.block[db]), 2);
1204                     fr.block[db].id=enxDISRE;
1205                     fr.block[db].sub[0].nr=ndisre;
1206                     fr.block[db].sub[1].nr=ndisre;
1207 #ifndef GMX_DOUBLE
1208                     fr.block[db].sub[0].type=xdr_datatype_float;
1209                     fr.block[db].sub[1].type=xdr_datatype_float;
1210                     fr.block[db].sub[0].fval=disre_rt;
1211                     fr.block[db].sub[1].fval=disre_rm3tav;
1212 #else
1213                     fr.block[db].sub[0].type=xdr_datatype_double;
1214                     fr.block[db].sub[1].type=xdr_datatype_double;
1215                     fr.block[db].sub[0].dval=disre_rt;
1216                     fr.block[db].sub[1].dval=disre_rm3tav;
1217 #endif
1218                 }
1219                 /* here we can put new-style blocks */
1220
1221                 /* Free energy perturbation blocks */
1222                 if (md->dhc)
1223                 {
1224                     mde_delta_h_coll_handle_block(md->dhc, &fr, fr.nblock);
1225                 }
1226
1227                 /* we can now free & reset the data in the blocks */
1228                 if (md->dhc)
1229                 {
1230                     mde_delta_h_coll_reset(md->dhc);
1231                 }
1232
1233                 /* do the actual I/O */
1234                 do_enx(fp_ene,&fr);
1235                 gmx_fio_check_file_position(enx_file_pointer(fp_ene));
1236                 if (fr.nre)
1237                 {
1238                     /* We have stored the sums, so reset the sum history */
1239                     reset_ebin_sums(md->ebin);
1240                 }
1241             }
1242             free_enxframe(&fr);
1243             break;
1244         case eprAVER:
1245             if (log)
1246             {
1247                 pprint(log,"A V E R A G E S",md);
1248             }
1249             break;
1250         case eprRMS:
1251             if (log)
1252             {
1253                 pprint(log,"R M S - F L U C T U A T I O N S",md);
1254             }
1255             break;
1256         default:
1257             gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
1258     }
1259
1260     if (log)
1261     {
1262         for(i=0;i<opts->ngtc;i++)
1263         {
1264             if(opts->annealing[i]!=eannNO)
1265             {
1266                 fprintf(log,"Current ref_t for group %s: %8.1f\n",
1267                         *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
1268                         opts->ref_t[i]);
1269             }
1270         }
1271         if (mode==eprNORMAL && fcd->orires.nr>0)
1272         {
1273             print_orires_log(log,&(fcd->orires));
1274         }
1275         fprintf(log,"   Energies (%s)\n",unit_energy);
1276         pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
1277         fprintf(log,"\n");
1278
1279         if (!bCompact)
1280         {
1281             if (md->bDynBox)
1282             {
1283                 pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
1284                         mode,TRUE);
1285                 fprintf(log,"\n");
1286             }
1287             if (md->bConstrVir)
1288             {
1289                 fprintf(log,"   Constraint Virial (%s)\n",unit_energy);
1290                 pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
1291                 fprintf(log,"\n");
1292                 fprintf(log,"   Force Virial (%s)\n",unit_energy);
1293                 pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
1294                 fprintf(log,"\n");
1295             }
1296             if (md->bVir)
1297             {
1298                 fprintf(log,"   Total Virial (%s)\n",unit_energy);
1299                 pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
1300                 fprintf(log,"\n");
1301             }
1302             if (md->bPress)
1303             {
1304                 fprintf(log,"   Pressure (%s)\n",unit_pres_bar);
1305                 pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
1306                 fprintf(log,"\n");
1307             }
1308             fprintf(log,"   Total Dipole (%s)\n",unit_dipole_D);
1309             pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
1310             fprintf(log,"\n");
1311
1312             if (md->nE > 1)
1313             {
1314                 if (md->print_grpnms==NULL)
1315                 {
1316                     snew(md->print_grpnms,md->nE);
1317                     n=0;
1318                     for(i=0; (i<md->nEg); i++)
1319                     {
1320                         ni=groups->grps[egcENER].nm_ind[i];
1321                         for(j=i; (j<md->nEg); j++)
1322                         {
1323                             nj=groups->grps[egcENER].nm_ind[j];
1324                             sprintf(buf,"%s-%s",*(groups->grpname[ni]),
1325                                     *(groups->grpname[nj]));
1326                             md->print_grpnms[n++]=strdup(buf);
1327                         }
1328                     }
1329                 }
1330                 sprintf(buf,"Epot (%s)",unit_energy);
1331                 fprintf(log,"%15s   ",buf);
1332                 for(i=0; (i<egNR); i++)
1333                 {
1334                     if (md->bEInd[i])
1335                     {
1336                         fprintf(log,"%12s   ",egrp_nm[i]);
1337                     }
1338                 }
1339                 fprintf(log,"\n");
1340                 for(i=0; (i<md->nE); i++)
1341                 {
1342                     fprintf(log,"%15s",md->print_grpnms[i]);
1343                     pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
1344                             FALSE);
1345                 }
1346                 fprintf(log,"\n");
1347             }
1348             if (md->nTC > 1)
1349             {
1350                 pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
1351                 fprintf(log,"\n");
1352             }
1353             if (md->nU > 1)
1354             {
1355                 fprintf(log,"%15s   %12s   %12s   %12s\n",
1356                         "Group","Ux","Uy","Uz");
1357                 for(i=0; (i<md->nU); i++)
1358                 {
1359                     ni=groups->grps[egcACC].nm_ind[i];
1360                     fprintf(log,"%15s",*groups->grpname[ni]);
1361                     pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
1362                 }
1363                 fprintf(log,"\n");
1364             }
1365         }
1366     }
1367
1368 }
1369
1370 void update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
1371 {
1372     int i;
1373
1374     enerhist->nsteps     = mdebin->ebin->nsteps;
1375     enerhist->nsum       = mdebin->ebin->nsum;
1376     enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
1377     enerhist->nsum_sim   = mdebin->ebin->nsum_sim;
1378     enerhist->nener      = mdebin->ebin->nener;
1379
1380     if (mdebin->ebin->nsum > 0)
1381     {
1382         /* Check if we need to allocate first */
1383         if(enerhist->ener_ave == NULL)
1384         {
1385             snew(enerhist->ener_ave,enerhist->nener);
1386             snew(enerhist->ener_sum,enerhist->nener);
1387         }
1388
1389         for(i=0;i<enerhist->nener;i++)
1390         {
1391             enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
1392             enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
1393         }
1394     }
1395
1396     if (mdebin->ebin->nsum_sim > 0)
1397     {
1398         /* Check if we need to allocate first */
1399         if(enerhist->ener_sum_sim == NULL)
1400         {
1401             snew(enerhist->ener_sum_sim,enerhist->nener);
1402         }
1403
1404         for(i=0;i<enerhist->nener;i++)
1405         {
1406             enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
1407         }
1408     }
1409     if (mdebin->dhc)
1410     {
1411         mde_delta_h_coll_update_energyhistory(mdebin->dhc, enerhist);
1412     }
1413 }
1414
1415 void restore_energyhistory_from_state(t_mdebin * mdebin,
1416                                       energyhistory_t * enerhist)
1417 {
1418     int i;
1419
1420     if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
1421         mdebin->ebin->nener != enerhist->nener)
1422     {
1423         gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1424                   mdebin->ebin->nener,enerhist->nener);
1425     }
1426
1427     mdebin->ebin->nsteps     = enerhist->nsteps;
1428     mdebin->ebin->nsum       = enerhist->nsum;
1429     mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
1430     mdebin->ebin->nsum_sim   = enerhist->nsum_sim;
1431
1432     for(i=0; i<mdebin->ebin->nener; i++)
1433     {
1434         mdebin->ebin->e[i].eav  =
1435                   (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
1436         mdebin->ebin->e[i].esum =
1437                   (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
1438         mdebin->ebin->e_sim[i].esum =
1439                   (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
1440     }
1441     if (mdebin->dhc)
1442     {
1443         mde_delta_h_coll_restore_energyhistory(mdebin->dhc, enerhist);
1444     }
1445 }