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