src/gromacs/mdlib/mdebin.c

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