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
  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
  88 t_mdebin *init_mdebin(ener_file_t fp_ene,
  89                       const gmx_mtop_t *mtop,
  90                       const t_inputrec *ir,
  91                       FILE *fp_dhdl)
  92 {
  93     const char *ener_nm[F_NRE];
  94     static const char *vir_nm[] = {
  95         "Vir-XX", "Vir-XY", "Vir-XZ",
  96         "Vir-YX", "Vir-YY", "Vir-YZ",
  97         "Vir-ZX", "Vir-ZY", "Vir-ZZ"
  98     };
  99     static const char *sv_nm[] = {
 100         "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
 101         "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
 102         "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
 103     };
 104     static const char *fv_nm[] = {
 105         "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
 106         "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
 107         "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
 108     };
 109     static const char *pres_nm[] = {
 110         "Pres-XX","Pres-XY","Pres-XZ",
 111         "Pres-YX","Pres-YY","Pres-YZ",
 112         "Pres-ZX","Pres-ZY","Pres-ZZ"
 113     };
 114     static const char *surft_nm[] = {
 115         "#Surf*SurfTen"
 116     };
 117     static const char *mu_nm[] = {
 118         "Mu-X", "Mu-Y", "Mu-Z"
 119     };
 120     static const char *vcos_nm[] = {
 121         "2CosZ*Vel-X"
 122     };
 123     static const char *visc_nm[] = {
 124         "1/Viscosity"
 125     };
 126     static const char *baro_nm[] = {
 127         "Barostat"
 128     };
 129
 130     char     **grpnms;
 131     const gmx_groups_t *groups;
 132     char     **gnm;
 133     char     buf[256];
 134     const char     *bufi;
 135     t_mdebin *md;
 136     int      i,j,ni,nj,n,nh,k,kk,ncon,nset;
 137     gmx_bool     bBHAM,bNoseHoover,b14;
 138
 139     snew(md,1);
 140
 141     if (EI_DYNAMICS(ir->eI))
 142     {
 143         md->delta_t = ir->delta_t;
 144     }
 145     else
 146     {
 147         md->delta_t = 0;
 148     }
 149
 150     groups = &mtop->groups;
 151
 152     bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
 153     b14   = (gmx_mtop_ftype_count(mtop,F_LJ14) > 0 ||
 154              gmx_mtop_ftype_count(mtop,F_LJC14_Q) > 0);
 155
 156     ncon = gmx_mtop_ftype_count(mtop,F_CONSTR);
 157     nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
 158     md->bConstr    = (ncon > 0 || nset > 0);
 159     md->bConstrVir = FALSE;
 160     if (md->bConstr) {
 161         if (ncon > 0 && ir->eConstrAlg == econtLINCS) {
 162             if (ir->eI == eiSD2)
 163                 md->nCrmsd = 2;
 164             else
 165                 md->nCrmsd = 1;
 166         }
 167         md->bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL);
 168     } else {
 169         md->nCrmsd = 0;
 170     }
 171
 172     /* Energy monitoring */
 173     for(i=0;i<egNR;i++)
 174     {
 175         md->bEInd[i]=FALSE;
 176     }
 177
 178 #ifndef GMX_OPENMM
 179     for(i=0; i<F_NRE; i++)
 180     {
 181         md->bEner[i] = FALSE;
 182         if (i == F_LJ)
 183             md->bEner[i] = !bBHAM;
 184         else if (i == F_BHAM)
 185             md->bEner[i] = bBHAM;
 186         else if (i == F_EQM)
 187             md->bEner[i] = ir->bQMMM;
 188         else if (i == F_COUL_LR)
 189             md->bEner[i] = (ir->rcoulomb > ir->rlist);
 190         else if (i == F_LJ_LR)
 191             md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist);
 192         else if (i == F_BHAM_LR)
 193             md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist);
 194         else if (i == F_RF_EXCL)
 195             md->bEner[i] = (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC);
 196         else if (i == F_COUL_RECIP)
 197             md->bEner[i] = EEL_FULL(ir->coulombtype);
 198         else if (i == F_LJ14)
 199             md->bEner[i] = b14;
 200         else if (i == F_COUL14)
 201             md->bEner[i] = b14;
 202         else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
 203             md->bEner[i] = FALSE;
 204         else if ((i == F_DVDL) || (i == F_DKDL))
 205             md->bEner[i] = (ir->efep != efepNO);
 206         else if (i == F_DHDL_CON)
 207             md->bEner[i] = (ir->efep != efepNO && md->bConstr);
 208         else if ((interaction_function[i].flags & IF_VSITE) ||
 209                  (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
 210             md->bEner[i] = FALSE;
 211         else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES)  || (i==F_EQM))
 212             md->bEner[i] = TRUE;
 213         else if ((i == F_GBPOL) && ir->implicit_solvent==eisGBSA)
 214             md->bEner[i] = TRUE;
 215         else if ((i == F_NPSOLVATION) && ir->implicit_solvent==eisGBSA && (ir->sa_algorithm != esaNO))
 216             md->bEner[i] = TRUE;
 217         else if ((i == F_GB12) || (i == F_GB13) || (i == F_GB14))
 218             md->bEner[i] = FALSE;
 219         else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
 220             md->bEner[i] = EI_DYNAMICS(ir->eI);
 221         else if (i==F_VTEMP)
 222             md->bEner[i] =  (EI_DYNAMICS(ir->eI) && getenv("GMX_VIRIAL_TEMPERATURE"));
 223         else if (i == F_DISPCORR || i == F_PDISPCORR)
 224             md->bEner[i] = (ir->eDispCorr != edispcNO);
 225         else if (i == F_DISRESVIOL)
 226             md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_DISRES) > 0);
 227         else if (i == F_ORIRESDEV)
 228             md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0);
 229         else if (i == F_CONNBONDS)
 230             md->bEner[i] = FALSE;
 231         else if (i == F_COM_PULL)
 232             md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F);
 233         else if (i == F_ECONSERVED)
 234             md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) &&
 235                             (ir->epc == epcNO || ir->epc==epcMTTK));
 236         else
 237             md->bEner[i] = (gmx_mtop_ftype_count(mtop,i) > 0);
 238     }
 239 #else
 240     /* OpenMM always produces only the following 4 energy terms */
 241     md->bEner[F_EPOT] = TRUE;
 242     md->bEner[F_EKIN] = TRUE;
 243     md->bEner[F_ETOT] = TRUE;
 244     md->bEner[F_TEMP] = TRUE;
 245 #endif
 246
 247     md->f_nre=0;
 248     for(i=0; i<F_NRE; i++)
 249     {
 250         if (md->bEner[i])
 251         {
 252             /* FIXME: The constness should not be cast away */
 253             /*ener_nm[f_nre]=(char *)interaction_function[i].longname;*/
 254             ener_nm[md->f_nre]=interaction_function[i].longname;
 255             md->f_nre++;
 256         }
 257     }
 258
 259     md->epc = ir->epc;
 260     md->bDiagPres = !TRICLINIC(ir->ref_p);
 261     md->ref_p = (ir->ref_p[XX][XX]+ir->ref_p[YY][YY]+ir->ref_p[ZZ][ZZ])/DIM;
 262     md->bTricl = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
 263     md->bDynBox = DYNAMIC_BOX(*ir);
 264     md->etc = ir->etc;
 265     md->bNHC_trotter = IR_NVT_TROTTER(ir);
 266     md->bMTTK = IR_NPT_TROTTER(ir);
 267
 268     md->ebin  = mk_ebin();
 269     /* Pass NULL for unit to let get_ebin_space determine the units
 270      * for interaction_function[i].longname
 271      */
 272     md->ie    = get_ebin_space(md->ebin,md->f_nre,ener_nm,NULL);
 273     if (md->nCrmsd)
 274     {
 275         /* This should be called directly after the call for md->ie,
 276          * such that md->iconrmsd follows directly in the list.
 277          */
 278         md->iconrmsd = get_ebin_space(md->ebin,md->nCrmsd,conrmsd_nm,"");
 279     }
 280     if (md->bDynBox)
 281     {
 282         md->ib    = get_ebin_space(md->ebin,
 283                                    md->bTricl ? NTRICLBOXS : NBOXS,
 284                                    md->bTricl ? tricl_boxs_nm : boxs_nm,
 285                                    unit_length);
 286         md->ivol  = get_ebin_space(md->ebin, 1, vol_nm,  unit_volume);
 287         md->idens = get_ebin_space(md->ebin, 1, dens_nm, unit_density_SI);
 288         if (md->bDiagPres)
 289         {
 290             md->ipv   = get_ebin_space(md->ebin, 1, pv_nm,   unit_energy);
 291             md->ienthalpy = get_ebin_space(md->ebin, 1, enthalpy_nm,   unit_energy);
 292         }
 293     }
 294     if (md->bConstrVir)
 295     {
 296         md->isvir = get_ebin_space(md->ebin,asize(sv_nm),sv_nm,unit_energy);
 297         md->ifvir = get_ebin_space(md->ebin,asize(fv_nm),fv_nm,unit_energy);
 298     }
 299     md->ivir   = get_ebin_space(md->ebin,asize(vir_nm),vir_nm,unit_energy);
 300     md->ipres  = get_ebin_space(md->ebin,asize(pres_nm),pres_nm,unit_pres_bar);
 301     md->isurft = get_ebin_space(md->ebin,asize(surft_nm),surft_nm,
 302                                 unit_surft_bar);
 303     if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
 304     {
 305         md->ipc = get_ebin_space(md->ebin,md->bTricl ? 6 : 3,
 306                                  boxvel_nm,unit_vel);
 307     }
 308     md->imu    = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
 309     if (ir->cos_accel != 0)
 310     {
 311         md->ivcos = get_ebin_space(md->ebin,asize(vcos_nm),vcos_nm,unit_vel);
 312         md->ivisc = get_ebin_space(md->ebin,asize(visc_nm),visc_nm,
 313                                    unit_invvisc_SI);
 314     }
 315
 316     /* Energy monitoring */
 317     for(i=0;i<egNR;i++)
 318     {
 319         md->bEInd[i] = FALSE;
 320     }
 321     md->bEInd[egCOULSR] = TRUE;
 322     md->bEInd[egLJSR  ] = TRUE;
 323
 324     if (ir->rcoulomb > ir->rlist)
 325     {
 326         md->bEInd[egCOULLR] = TRUE;
 327     }
 328     if (!bBHAM)
 329     {
 330         if (ir->rvdw > ir->rlist)
 331         {
 332             md->bEInd[egLJLR]   = TRUE;
 333         }
 334     }
 335     else
 336     {
 337         md->bEInd[egLJSR]   = FALSE;
 338         md->bEInd[egBHAMSR] = TRUE;
 339         if (ir->rvdw > ir->rlist)
 340         {
 341             md->bEInd[egBHAMLR]   = TRUE;
 342         }
 343     }
 344     if (b14)
 345     {
 346         md->bEInd[egLJ14] = TRUE;
 347         md->bEInd[egCOUL14] = TRUE;
 348     }
 349     md->nEc=0;
 350     for(i=0; (i<egNR); i++)
 351     {
 352         if (md->bEInd[i])
 353         {
 354             md->nEc++;
 355         }
 356     }
 357
 358     n=groups->grps[egcENER].nr;
 359     md->nEg=n;
 360     md->nE=(n*(n+1))/2;
 361     snew(md->igrp,md->nE);
 362     if (md->nE > 1)
 363     {
 364         n=0;
 365         snew(gnm,md->nEc);
 366         for(k=0; (k<md->nEc); k++)
 367         {
 368             snew(gnm[k],STRLEN);
 369         }
 370         for(i=0; (i<groups->grps[egcENER].nr); i++)
 371         {
 372             ni=groups->grps[egcENER].nm_ind[i];
 373             for(j=i; (j<groups->grps[egcENER].nr); j++)
 374             {
 375                 nj=groups->grps[egcENER].nm_ind[j];
 376                 for(k=kk=0; (k<egNR); k++)
 377                 {
 378                     if (md->bEInd[k])
 379                     {
 380                         sprintf(gnm[kk],"%s:%s-%s",egrp_nm[k],
 381                                 *(groups->grpname[ni]),*(groups->grpname[nj]));
 382                         kk++;
 383                     }
 384                 }
 385                 md->igrp[n]=get_ebin_space(md->ebin,md->nEc,
 386                                            (const char **)gnm,unit_energy);
 387                 n++;
 388             }
 389         }
 390         for(k=0; (k<md->nEc); k++)
 391         {
 392             sfree(gnm[k]);
 393         }
 394         sfree(gnm);
 395
 396         if (n != md->nE)
 397         {
 398             gmx_incons("Number of energy terms wrong");
 399         }
 400     }
 401
 402     md->nTC=groups->grps[egcTC].nr;
 403     md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */
 404     if (md->bMTTK)
 405     {
 406         md->nTCP = 1;  /* assume only one possible coupling system for barostat
 407                           for now */
 408     }
 409     else
 410     {
 411         md->nTCP = 0;
 412     }
 413
 414     if (md->etc == etcNOSEHOOVER) {
 415         if (md->bNHC_trotter) {
 416             md->mde_n = 2*md->nNHC*md->nTC;
 417         }
 418         else
 419         {
 420             md->mde_n = 2*md->nTC;
 421         }
 422         if (md->epc == epcMTTK)
 423         {
 424             md->mdeb_n = 2*md->nNHC*md->nTCP;
 425         }
 426     } else {
 427         md->mde_n = md->nTC;
 428         md->mdeb_n = 0;
 429     }
 430
 431     snew(md->tmp_r,md->mde_n);
 432     snew(md->tmp_v,md->mde_n);
 433     snew(md->grpnms,md->mde_n);
 434     grpnms = md->grpnms;
 435
 436     for(i=0; (i<md->nTC); i++)
 437     {
 438         ni=groups->grps[egcTC].nm_ind[i];
 439         sprintf(buf,"T-%s",*(groups->grpname[ni]));
 440         grpnms[i]=strdup(buf);
 441     }
 442     md->itemp=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,
 443                              unit_temp_K);
 444
 445     bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL); /* whether to print Nose-Hoover chains */
 446
 447     if (md->etc == etcNOSEHOOVER)
 448     {
 449         if (bNoseHoover)
 450         {
 451             if (md->bNHC_trotter)
 452             {
 453                 for(i=0; (i<md->nTC); i++)
 454                 {
 455                     ni=groups->grps[egcTC].nm_ind[i];
 456                     bufi = *(groups->grpname[ni]);
 457                     for(j=0; (j<md->nNHC); j++)
 458                     {
 459                         sprintf(buf,"Xi-%d-%s",j,bufi);
 460                         grpnms[2*(i*md->nNHC+j)]=strdup(buf);
 461                         sprintf(buf,"vXi-%d-%s",j,bufi);
 462                         grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
 463                     }
 464                 }
 465                 md->itc=get_ebin_space(md->ebin,md->mde_n,
 466                                        (const char **)grpnms,unit_invtime);
 467                 if (md->bMTTK)
 468                 {
 469                     for(i=0; (i<md->nTCP); i++)
 470                     {
 471                         bufi = baro_nm[0];  /* All barostat DOF's together for now. */
 472                         for(j=0; (j<md->nNHC); j++)
 473                         {
 474                             sprintf(buf,"Xi-%d-%s",j,bufi);
 475                             grpnms[2*(i*md->nNHC+j)]=strdup(buf);
 476                             sprintf(buf,"vXi-%d-%s",j,bufi);
 477                             grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
 478                         }
 479                     }
 480                     md->itcb=get_ebin_space(md->ebin,md->mdeb_n,
 481                                             (const char **)grpnms,unit_invtime);
 482                 }
 483             }
 484             else
 485             {
 486                 for(i=0; (i<md->nTC); i++)
 487                 {
 488                     ni=groups->grps[egcTC].nm_ind[i];
 489                     bufi = *(groups->grpname[ni]);
 490                     sprintf(buf,"Xi-%s",bufi);
 491                     grpnms[2*i]=strdup(buf);
 492                     sprintf(buf,"vXi-%s",bufi);
 493                     grpnms[2*i+1]=strdup(buf);
 494                 }
 495                 md->itc=get_ebin_space(md->ebin,md->mde_n,
 496                                        (const char **)grpnms,unit_invtime);
 497             }
 498         }
 499     }
 500     else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
 501              md->etc == etcVRESCALE)
 502     {
 503         for(i=0; (i<md->nTC); i++)
 504         {
 505             ni=groups->grps[egcTC].nm_ind[i];
 506             sprintf(buf,"Lamb-%s",*(groups->grpname[ni]));
 507             grpnms[i]=strdup(buf);
 508         }
 509         md->itc=get_ebin_space(md->ebin,md->mde_n,(const char **)grpnms,"");
 510     }
 511
 512     sfree(grpnms);
 513
 514
 515     md->nU=groups->grps[egcACC].nr;
 516     if (md->nU > 1)
 517     {
 518         snew(grpnms,3*md->nU);
 519         for(i=0; (i<md->nU); i++)
 520         {
 521             ni=groups->grps[egcACC].nm_ind[i];
 522             sprintf(buf,"Ux-%s",*(groups->grpname[ni]));
 523             grpnms[3*i+XX]=strdup(buf);
 524             sprintf(buf,"Uy-%s",*(groups->grpname[ni]));
 525             grpnms[3*i+YY]=strdup(buf);
 526             sprintf(buf,"Uz-%s",*(groups->grpname[ni]));
 527             grpnms[3*i+ZZ]=strdup(buf);
 528         }
 529         md->iu=get_ebin_space(md->ebin,3*md->nU,(const char **)grpnms,unit_vel);
 530         sfree(grpnms);
 531     }
 532
 533     if ( fp_ene )
 534     {
 535         do_enxnms(fp_ene,&md->ebin->nener,&md->ebin->enm);
 536     }
 537
 538     md->print_grpnms=NULL;
 539
 540     /* check whether we're going to write dh histograms */
 541     md->dhc=NULL;
 542     if (ir->separate_dhdl_file == sepdhdlfileNO )
 543     {
 544         int i;
 545         snew(md->dhc, 1);
 546
 547         mde_delta_h_coll_init(md->dhc, ir);
 548         md->fp_dhdl = NULL;
 549     }
 550     else
 551     {
 552         md->fp_dhdl = fp_dhdl;
 553     }
 554     md->dhdl_derivatives = (ir->dhdl_derivatives==dhdlderivativesYES);
 555     return md;
 556 }
 557
 558 FILE *open_dhdl(const char *filename,const t_inputrec *ir,
 559                 const output_env_t oenv)
 560 {
 561     FILE *fp;
 562     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
 563     char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
 564     char **setname;
 565     char buf[STRLEN];
 566
 567     sprintf(label_x,"%s (%s)","Time",unit_time);
 568     if (ir->n_flambda == 0)
 569     {
 570         sprintf(title,"%s",dhdl);
 571         sprintf(label_y,"%s (%s %s)",
 572                 dhdl,unit_energy,"[\\lambda]\\S-1\\N");
 573     }
 574     else
 575     {
 576         sprintf(title,"%s, %s",dhdl,deltag);
 577         sprintf(label_y,"(%s)",unit_energy);
 578     }
 579     fp = gmx_fio_fopen(filename,"w+");
 580     xvgr_header(fp,title,label_x,label_y,exvggtXNY,oenv);
 581
 582     if (ir->delta_lambda == 0)
 583     {
 584         sprintf(buf,"T = %g (K), %s = %g",
 585                 ir->opts.ref_t[0],lambda,ir->init_lambda);
 586     }
 587     else
 588     {
 589         sprintf(buf,"T = %g (K)",
 590                 ir->opts.ref_t[0]);
 591     }
 592     xvgr_subtitle(fp,buf,oenv);
 593
 594     if (ir->n_flambda > 0)
 595     {
 596         int nsets,s,nsi=0;
 597         /* g_bar has to determine the lambda values used in this simulation
 598          * from this xvg legend.  */
 599         nsets = ( (ir->dhdl_derivatives==dhdlderivativesYES) ? 1 : 0) +
 600                   ir->n_flambda;
 601         snew(setname,nsets);
 602         if (ir->dhdl_derivatives == dhdlderivativesYES)
 603         {
 604             sprintf(buf,"%s %s %g",dhdl,lambda,ir->init_lambda);
 605             setname[nsi++] = gmx_strdup(buf);
 606         }
 607         for(s=0; s<ir->n_flambda; s++)
 608         {
 609             sprintf(buf,"%s %s %g",deltag,lambda,ir->flambda[s]);
 610             setname[nsi++] = gmx_strdup(buf);
 611         }
 612         xvgr_legend(fp,nsets,(const char**)setname,oenv);
 613
 614         for(s=0; s<nsets; s++)
 615         {
 616             sfree(setname[s]);
 617         }
 618         sfree(setname);
 619     }
 620
 621     return fp;
 622 }
 623
 624 static void copy_energy(t_mdebin *md, real e[],real ecpy[])
 625 {
 626     int i,j;
 627
 628     for(i=j=0; (i<F_NRE); i++)
 629         if (md->bEner[i])
 630             ecpy[j++] = e[i];
 631     if (j != md->f_nre)
 632         gmx_incons("Number of energy terms wrong");
 633 }
 634
 635 void upd_mdebin(t_mdebin *md, gmx_bool write_dhdl,
 636                 gmx_bool bSum,
 637                 double time,
 638                 real tmass,
 639                 gmx_enerdata_t *enerd,
 640                 t_state *state,
 641                 matrix  box,
 642                 tensor svir,
 643                 tensor fvir,
 644                 tensor vir,
 645                 tensor pres,
 646                 gmx_ekindata_t *ekind,
 647                 rvec mu_tot,
 648                 gmx_constr_t constr)
 649 {
 650     int    i,j,k,kk,m,n,gid;
 651     real   crmsd[2],tmp6[6];
 652     real   bs[NTRICLBOXS],vol,dens,pv,enthalpy;
 653     real   eee[egNR];
 654     real   ecopy[F_NRE];
 655     real   tmp;
 656     gmx_bool   bNoseHoover;
 657
 658     /* Do NOT use the box in the state variable, but the separate box provided
 659      * as an argument. This is because we sometimes need to write the box from
 660      * the last timestep to match the trajectory frames.
 661      */
 662     copy_energy(md, enerd->term,ecopy);
 663     add_ebin(md->ebin,md->ie,md->f_nre,ecopy,bSum);
 664     if (md->nCrmsd)
 665     {
 666         crmsd[0] = constr_rmsd(constr,FALSE);
 667         if (md->nCrmsd > 1)
 668         {
 669             crmsd[1] = constr_rmsd(constr,TRUE);
 670         }
 671         add_ebin(md->ebin,md->iconrmsd,md->nCrmsd,crmsd,FALSE);
 672     }
 673     if (md->bDynBox)
 674     {
 675         int nboxs;
 676         if(md->bTricl)
 677         {
 678             bs[0] = box[XX][XX];
 679             bs[1] = box[YY][YY];
 680             bs[2] = box[ZZ][ZZ];
 681             bs[3] = box[YY][XX];
 682             bs[4] = box[ZZ][XX];
 683             bs[5] = box[ZZ][YY];
 684             nboxs=NTRICLBOXS;
 685         }
 686         else
 687         {
 688             bs[0] = box[XX][XX];
 689             bs[1] = box[YY][YY];
 690             bs[2] = box[ZZ][ZZ];
 691             nboxs=NBOXS;
 692         }
 693         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 694         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 695
 696         add_ebin(md->ebin,md->ib   ,nboxs,bs   ,bSum);
 697         add_ebin(md->ebin,md->ivol ,1    ,&vol ,bSum);
 698         add_ebin(md->ebin,md->idens,1    ,&dens,bSum);
 699
 700         if (md->bDiagPres)
 701         {
 702             /* This is pV (in kJ/mol).  The pressure is the reference pressure,
 703                not the instantaneous pressure */
 704             pv = vol*md->ref_p/PRESFAC;
 705
 706             add_ebin(md->ebin,md->ipv  ,1    ,&pv  ,bSum);
 707             enthalpy = pv + enerd->term[F_ETOT];
 708             add_ebin(md->ebin,md->ienthalpy  ,1    ,&enthalpy  ,bSum);
 709         }
 710     }
 711     if (md->bConstrVir)
 712     {
 713         add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
 714         add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
 715     }
 716     add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
 717     add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
 718     tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
 719     add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
 720     if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
 721     {
 722         tmp6[0] = state->boxv[XX][XX];
 723         tmp6[1] = state->boxv[YY][YY];
 724         tmp6[2] = state->boxv[ZZ][ZZ];
 725         tmp6[3] = state->boxv[YY][XX];
 726         tmp6[4] = state->boxv[ZZ][XX];
 727         tmp6[5] = state->boxv[ZZ][YY];
 728         add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
 729     }
 730     add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
 731     if (ekind && ekind->cosacc.cos_accel != 0)
 732     {
 733         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 734         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 735         add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
 736         /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
 737         tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
 738                  *dens*vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
 739         add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
 740     }
 741     if (md->nE > 1)
 742     {
 743         n=0;
 744         for(i=0; (i<md->nEg); i++)
 745         {
 746             for(j=i; (j<md->nEg); j++)
 747             {
 748                 gid=GID(i,j,md->nEg);
 749                 for(k=kk=0; (k<egNR); k++)
 750                 {
 751                     if (md->bEInd[k])
 752                     {
 753                         eee[kk++] = enerd->grpp.ener[k][gid];
 754                     }
 755                 }
 756                 add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
 757                 n++;
 758             }
 759         }
 760     }
 761
 762     if (ekind)
 763     {
 764         for(i=0; (i<md->nTC); i++)
 765         {
 766             md->tmp_r[i] = ekind->tcstat[i].T;
 767         }
 768         add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
 769
 770         /* whether to print Nose-Hoover chains: */
 771         bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL);
 772
 773         if (md->etc == etcNOSEHOOVER)
 774         {
 775             if (bNoseHoover)
 776             {
 777                 if (md->bNHC_trotter)
 778                 {
 779                     for(i=0; (i<md->nTC); i++)
 780                     {
 781                         for (j=0;j<md->nNHC;j++)
 782                         {
 783                             k = i*md->nNHC+j;
 784                             md->tmp_r[2*k] = state->nosehoover_xi[k];
 785                             md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
 786                         }
 787                     }
 788                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
 789
 790                     if (md->bMTTK) {
 791                         for(i=0; (i<md->nTCP); i++)
 792                         {
 793                             for (j=0;j<md->nNHC;j++)
 794                             {
 795                                 k = i*md->nNHC+j;
 796                                 md->tmp_r[2*k] = state->nhpres_xi[k];
 797                                 md->tmp_r[2*k+1] = state->nhpres_vxi[k];
 798                             }
 799                         }
 800                         add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
 801                     }
 802
 803                 }
 804                 else
 805                 {
 806                     for(i=0; (i<md->nTC); i++)
 807                     {
 808                         md->tmp_r[2*i] = state->nosehoover_xi[i];
 809                         md->tmp_r[2*i+1] = state->nosehoover_vxi[i];
 810                     }
 811                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
 812                 }
 813             }
 814         }
 815         else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
 816                  md->etc == etcVRESCALE)
 817         {
 818             for(i=0; (i<md->nTC); i++)
 819             {
 820                 md->tmp_r[i] = ekind->tcstat[i].lambda;
 821             }
 822             add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
 823         }
 824     }
 825
 826     if (ekind && md->nU > 1)
 827     {
 828         for(i=0; (i<md->nU); i++)
 829         {
 830             copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
 831         }
 832         add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
 833     }
 834
 835     ebin_increase_count(md->ebin,bSum);
 836
 837     /* BAR + thermodynamic integration values */
 838     if (write_dhdl)
 839     {
 840         if (md->fp_dhdl)
 841         {
 842             fprintf(md->fp_dhdl,"%.4f", time);
 843
 844             if (md->dhdl_derivatives)
 845             {
 846                 fprintf(md->fp_dhdl," %g", enerd->term[F_DVDL]+
 847                                            enerd->term[F_DKDL]+
 848                                            enerd->term[F_DHDL_CON]);
 849             }
 850             for(i=1; i<enerd->n_lambda; i++)
 851             {
 852                 fprintf(md->fp_dhdl," %g",
 853                         enerd->enerpart_lambda[i]-enerd->enerpart_lambda[0]);
 854             }
 855             fprintf(md->fp_dhdl,"\n");
 856         }
 857         /* and the binary BAR output */
 858         if (md->dhc)
 859         {
 860             mde_delta_h_coll_add_dh(md->dhc,
 861                                     enerd->term[F_DVDL]+ enerd->term[F_DKDL]+
 862                                     enerd->term[F_DHDL_CON],
 863                                     enerd->enerpart_lambda, time,
 864                                     state->lambda);
 865         }
 866     }
 867 }
 868
 869 void upd_mdebin_step(t_mdebin *md)
 870 {
 871     ebin_increase_count(md->ebin,FALSE);
 872 }
 873
 874 static void npr(FILE *log,int n,char c)
 875 {
 876     for(; (n>0); n--) fprintf(log,"%c",c);
 877 }
 878
 879 static void pprint(FILE *log,const char *s,t_mdebin *md)
 880 {
 881     char CHAR='#';
 882     int  slen;
 883     char buf1[22],buf2[22];
 884
 885     slen = strlen(s);
 886     fprintf(log,"\t<======  ");
 887     npr(log,slen,CHAR);
 888     fprintf(log,"  ==>\n");
 889     fprintf(log,"\t<====  %s  ====>\n",s);
 890     fprintf(log,"\t<==  ");
 891     npr(log,slen,CHAR);
 892     fprintf(log,"  ======>\n\n");
 893
 894     fprintf(log,"\tStatistics over %s steps using %s frames\n",
 895             gmx_step_str(md->ebin->nsteps_sim,buf1),
 896             gmx_step_str(md->ebin->nsum_sim,buf2));
 897     fprintf(log,"\n");
 898 }
 899
 900 void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lamb)
 901 {
 902     char buf[22];
 903
 904     fprintf(log,"   %12s   %12s   %12s\n"
 905             "   %12s   %12.5f   %12.5f\n\n",
 906             "Step","Time","Lambda",gmx_step_str(steps,buf),time,lamb);
 907 }
 908
 909 void print_ebin(ener_file_t fp_ene,gmx_bool bEne,gmx_bool bDR,gmx_bool bOR,
 910                 FILE *log,
 911                 gmx_large_int_t step,double time,
 912                 int mode,gmx_bool bCompact,
 913                 t_mdebin *md,t_fcdata *fcd,
 914                 gmx_groups_t *groups,t_grpopts *opts)
 915 {
 916     /*static char **grpnms=NULL;*/
 917     char        buf[246];
 918     int         i,j,n,ni,nj,ndr,nor,b;
 919     int         ndisre=0;
 920     real        *disre_rm3tav, *disre_rt;
 921
 922     /* these are for the old-style blocks (1 subblock, only reals), because
 923        there can be only one per ID for these */
 924     int         nr[enxNR];
 925     int         id[enxNR];
 926     real        *block[enxNR];
 927
 928     /* temporary arrays for the lambda values to write out */
 929     double      enxlambda_data[2];
 930
 931     t_enxframe  fr;
 932
 933     switch (mode)
 934     {
 935         case eprNORMAL:
 936             init_enxframe(&fr);
 937             fr.t            = time;
 938             fr.step         = step;
 939             fr.nsteps       = md->ebin->nsteps;
 940             fr.dt           = md->delta_t;
 941             fr.nsum         = md->ebin->nsum;
 942             fr.nre          = (bEne) ? md->ebin->nener : 0;
 943             fr.ener         = md->ebin->e;
 944             ndisre          = bDR ? fcd->disres.npair : 0;
 945             disre_rm3tav    = fcd->disres.rm3tav;
 946             disre_rt        = fcd->disres.rt;
 947             /* Optional additional old-style (real-only) blocks. */
 948             for(i=0; i<enxNR; i++)
 949             {
 950                 nr[i] = 0;
 951             }
 952             if (fcd->orires.nr > 0 && bOR)
 953             {
 954                 diagonalize_orires_tensors(&(fcd->orires));
 955                 nr[enxOR]     = fcd->orires.nr;
 956                 block[enxOR]  = fcd->orires.otav;
 957                 id[enxOR]     = enxOR;
 958                 nr[enxORI]    = (fcd->orires.oinsl != fcd->orires.otav) ?
 959                           fcd->orires.nr : 0;
 960                 block[enxORI] = fcd->orires.oinsl;
 961                 id[enxORI]    = enxORI;
 962                 nr[enxORT]    = fcd->orires.nex*12;
 963                 block[enxORT] = fcd->orires.eig;
 964                 id[enxORT]    = enxORT;
 965             }
 966
 967             /* whether we are going to wrte anything out: */
 968             if (fr.nre || ndisre || nr[enxOR] || nr[enxORI])
 969             {
 970
 971                 /* the old-style blocks go first */
 972                 fr.nblock = 0;
 973                 for(i=0; i<enxNR; i++)
 974                 {
 975                     if (nr[i] > 0)
 976                     {
 977                         fr.nblock = i + 1;
 978                     }
 979                 }
 980                 add_blocks_enxframe(&fr, fr.nblock);
 981                 for(b=0;b<fr.nblock;b++)
 982                 {
 983                     add_subblocks_enxblock(&(fr.block[b]), 1);
 984                     fr.block[b].id=id[b];
 985                     fr.block[b].sub[0].nr = nr[b];
 986 #ifndef GMX_DOUBLE
 987                     fr.block[b].sub[0].type = xdr_datatype_float;
 988                     fr.block[b].sub[0].fval = block[b];
 989 #else
 990                     fr.block[b].sub[0].type = xdr_datatype_double;
 991                     fr.block[b].sub[0].dval = block[b];
 992 #endif
 993                 }
 994
 995                 /* check for disre block & fill it. */
 996                 if (ndisre>0)
 997                 {
 998                     int db = fr.nblock;
 999                     fr.nblock+=1;
1000                     add_blocks_enxframe(&fr, fr.nblock);
1001
1002                     add_subblocks_enxblock(&(fr.block[db]), 2);
1003                     fr.block[db].id=enxDISRE;
1004                     fr.block[db].sub[0].nr=ndisre;
1005                     fr.block[db].sub[1].nr=ndisre;
1006 #ifndef GMX_DOUBLE
1007                     fr.block[db].sub[0].type=xdr_datatype_float;
1008                     fr.block[db].sub[1].type=xdr_datatype_float;
1009                     fr.block[db].sub[0].fval=disre_rt;
1010                     fr.block[db].sub[1].fval=disre_rm3tav;
1011 #else
1012                     fr.block[db].sub[0].type=xdr_datatype_double;
1013                     fr.block[db].sub[1].type=xdr_datatype_double;
1014                     fr.block[db].sub[0].dval=disre_rt;
1015                     fr.block[db].sub[1].dval=disre_rm3tav;
1016 #endif
1017                 }
1018                 /* here we can put new-style blocks */
1019
1020                 /* Free energy perturbation blocks */
1021                 if (md->dhc)
1022                 {
1023                     mde_delta_h_coll_handle_block(md->dhc, &fr, fr.nblock);
1024                 }
1025
1026                 /* do the actual I/O */
1027                 do_enx(fp_ene,&fr);
1028                 gmx_fio_check_file_position(enx_file_pointer(fp_ene));
1029                 if (fr.nre)
1030                 {
1031                     /* We have stored the sums, so reset the sum history */
1032                     reset_ebin_sums(md->ebin);
1033                 }
1034
1035                 /* we can now free & reset the data in the blocks */
1036                 if (md->dhc)
1037                     mde_delta_h_coll_reset(md->dhc);
1038             }
1039             free_enxframe(&fr);
1040             break;
1041         case eprAVER:
1042             if (log)
1043             {
1044                 pprint(log,"A V E R A G E S",md);
1045             }
1046             break;
1047         case eprRMS:
1048             if (log)
1049             {
1050                 pprint(log,"R M S - F L U C T U A T I O N S",md);
1051             }
1052             break;
1053         default:
1054             gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
1055     }
1056
1057     if (log)
1058     {
1059         for(i=0;i<opts->ngtc;i++)
1060         {
1061             if(opts->annealing[i]!=eannNO)
1062             {
1063                 fprintf(log,"Current ref_t for group %s: %8.1f\n",
1064                         *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
1065                         opts->ref_t[i]);
1066             }
1067         }
1068         if (mode==eprNORMAL && fcd->orires.nr>0)
1069         {
1070             print_orires_log(log,&(fcd->orires));
1071         }
1072         fprintf(log,"   Energies (%s)\n",unit_energy);
1073         pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
1074         fprintf(log,"\n");
1075
1076         if (!bCompact)
1077         {
1078             if (md->bDynBox)
1079             {
1080                 pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
1081                         mode,TRUE);
1082                 fprintf(log,"\n");
1083             }
1084             if (md->bConstrVir)
1085             {
1086                 fprintf(log,"   Constraint Virial (%s)\n",unit_energy);
1087                 pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
1088                 fprintf(log,"\n");
1089                 fprintf(log,"   Force Virial (%s)\n",unit_energy);
1090                 pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
1091                 fprintf(log,"\n");
1092             }
1093             fprintf(log,"   Total Virial (%s)\n",unit_energy);
1094             pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
1095             fprintf(log,"\n");
1096             fprintf(log,"   Pressure (%s)\n",unit_pres_bar);
1097             pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
1098             fprintf(log,"\n");
1099             fprintf(log,"   Total Dipole (%s)\n",unit_dipole_D);
1100             pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
1101             fprintf(log,"\n");
1102
1103             if (md->nE > 1)
1104             {
1105                 if (md->print_grpnms==NULL)
1106                 {
1107                     snew(md->print_grpnms,md->nE);
1108                     n=0;
1109                     for(i=0; (i<md->nEg); i++)
1110                     {
1111                         ni=groups->grps[egcENER].nm_ind[i];
1112                         for(j=i; (j<md->nEg); j++)
1113                         {
1114                             nj=groups->grps[egcENER].nm_ind[j];
1115                             sprintf(buf,"%s-%s",*(groups->grpname[ni]),
1116                                     *(groups->grpname[nj]));
1117                             md->print_grpnms[n++]=strdup(buf);
1118                         }
1119                     }
1120                 }
1121                 sprintf(buf,"Epot (%s)",unit_energy);
1122                 fprintf(log,"%15s   ",buf);
1123                 for(i=0; (i<egNR); i++)
1124                 {
1125                     if (md->bEInd[i])
1126                     {
1127                         fprintf(log,"%12s   ",egrp_nm[i]);
1128                     }
1129                 }
1130                 fprintf(log,"\n");
1131                 for(i=0; (i<md->nE); i++)
1132                 {
1133                     fprintf(log,"%15s",md->print_grpnms[i]);
1134                     pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
1135                             FALSE);
1136                 }
1137                 fprintf(log,"\n");
1138             }
1139             if (md->nTC > 1)
1140             {
1141                 pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
1142                 fprintf(log,"\n");
1143             }
1144             if (md->nU > 1)
1145             {
1146                 fprintf(log,"%15s   %12s   %12s   %12s\n",
1147                         "Group","Ux","Uy","Uz");
1148                 for(i=0; (i<md->nU); i++)
1149                 {
1150                     ni=groups->grps[egcACC].nm_ind[i];
1151                     fprintf(log,"%15s",*groups->grpname[ni]);
1152                     pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
1153                 }
1154                 fprintf(log,"\n");
1155             }
1156         }
1157     }
1158
1159 }
1160
1161 void update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
1162 {
1163     int i;
1164
1165     enerhist->nsteps     = mdebin->ebin->nsteps;
1166     enerhist->nsum       = mdebin->ebin->nsum;
1167     enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
1168     enerhist->nsum_sim   = mdebin->ebin->nsum_sim;
1169     enerhist->nener      = mdebin->ebin->nener;
1170
1171     if (mdebin->ebin->nsum > 0)
1172     {
1173         /* Check if we need to allocate first */
1174         if(enerhist->ener_ave == NULL)
1175         {
1176             snew(enerhist->ener_ave,enerhist->nener);
1177             snew(enerhist->ener_sum,enerhist->nener);
1178         }
1179
1180         for(i=0;i<enerhist->nener;i++)
1181         {
1182             enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
1183             enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
1184         }
1185     }
1186
1187     if (mdebin->ebin->nsum_sim > 0)
1188     {
1189         /* Check if we need to allocate first */
1190         if(enerhist->ener_sum_sim == NULL)
1191         {
1192             snew(enerhist->ener_sum_sim,enerhist->nener);
1193         }
1194
1195         for(i=0;i<enerhist->nener;i++)
1196         {
1197             enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
1198         }
1199     }
1200     if (mdebin->dhc)
1201     {
1202         mde_delta_h_coll_update_energyhistory(mdebin->dhc, enerhist);
1203     }
1204 }
1205
1206 void restore_energyhistory_from_state(t_mdebin * mdebin,
1207                                       energyhistory_t * enerhist)
1208 {
1209     int i;
1210
1211     if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
1212         mdebin->ebin->nener != enerhist->nener)
1213     {
1214         gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1215                   mdebin->ebin->nener,enerhist->nener);
1216     }
1217
1218     mdebin->ebin->nsteps     = enerhist->nsteps;
1219     mdebin->ebin->nsum       = enerhist->nsum;
1220     mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
1221     mdebin->ebin->nsum_sim   = enerhist->nsum_sim;
1222
1223     for(i=0; i<mdebin->ebin->nener; i++)
1224     {
1225         mdebin->ebin->e[i].eav  =
1226                   (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
1227         mdebin->ebin->e[i].esum =
1228                   (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
1229         mdebin->ebin->e_sim[i].esum =
1230                   (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
1231     }
1232     if (mdebin->dhc)
1233     {
1234         mde_delta_h_coll_restore_energyhistory(mdebin->dhc, enerhist);
1235     }
1236 }