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