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