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     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     md->nTC=groups->grps[egcTC].nr;
 399     md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */
 400     if (md->bMTTK)
 401     {
 402         md->nTCP = 1;  /* assume only one possible coupling system for barostat
 403                           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,
 462                                        (const char **)grpnms,unit_invtime);
 463                 if (md->bMTTK)
 464                 {
 465                     for(i=0; (i<md->nTCP); i++)
 466                     {
 467                         bufi = baro_nm[0];  /* All barostat DOF's together for now. */
 468                         for(j=0; (j<md->nNHC); j++)
 469                         {
 470                             sprintf(buf,"Xi-%d-%s",j,bufi);
 471                             grpnms[2*(i*md->nNHC+j)]=strdup(buf);
 472                             sprintf(buf,"vXi-%d-%s",j,bufi);
 473                             grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
 474                         }
 475                     }
 476                     md->itcb=get_ebin_space(md->ebin,md->mdeb_n,
 477                                             (const char **)grpnms,unit_invtime);
 478                 }
 479             }
 480             else
 481             {
 482                 for(i=0; (i<md->nTC); i++)
 483                 {
 484                     ni=groups->grps[egcTC].nm_ind[i];
 485                     bufi = *(groups->grpname[ni]);
 486                     sprintf(buf,"Xi-%s",bufi);
 487                     grpnms[2*i]=strdup(buf);
 488                     sprintf(buf,"vXi-%s",bufi);
 489                     grpnms[2*i+1]=strdup(buf);
 490                 }
 491                 md->itc=get_ebin_space(md->ebin,md->mde_n,
 492                                        (const char **)grpnms,unit_invtime);
 493             }
 494         }
 495     }
 496     else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
 497              md->etc == etcVRESCALE)
 498     {
 499         for(i=0; (i<md->nTC); i++)
 500         {
 501             ni=groups->grps[egcTC].nm_ind[i];
 502             sprintf(buf,"Lamb-%s",*(groups->grpname[ni]));
 503             grpnms[i]=strdup(buf);
 504         }
 505         md->itc=get_ebin_space(md->ebin,md->mde_n,(const char **)grpnms,"");
 506     }
 507
 508     sfree(grpnms);
 509
 510
 511     md->nU=groups->grps[egcACC].nr;
 512     if (md->nU > 1)
 513     {
 514         snew(grpnms,3*md->nU);
 515         for(i=0; (i<md->nU); i++)
 516         {
 517             ni=groups->grps[egcACC].nm_ind[i];
 518             sprintf(buf,"Ux-%s",*(groups->grpname[ni]));
 519             grpnms[3*i+XX]=strdup(buf);
 520             sprintf(buf,"Uy-%s",*(groups->grpname[ni]));
 521             grpnms[3*i+YY]=strdup(buf);
 522             sprintf(buf,"Uz-%s",*(groups->grpname[ni]));
 523             grpnms[3*i+ZZ]=strdup(buf);
 524         }
 525         md->iu=get_ebin_space(md->ebin,3*md->nU,(const char **)grpnms,unit_vel);
 526         sfree(grpnms);
 527     }
 528
 529     if ( fp_ene )
 530     {
 531         do_enxnms(fp_ene,&md->ebin->nener,&md->ebin->enm);
 532     }
 533
 534     md->print_grpnms=NULL;
 535
 536     /* check whether we're going to write dh histograms */
 537     md->dhc=NULL;
 538     if (ir->separate_dhdl_file == sepdhdlfileNO )
 539     {
 540         int i;
 541         snew(md->dhc, 1);
 542
 543         mde_delta_h_coll_init(md->dhc, ir);
 544         md->fp_dhdl = NULL;
 545     }
 546     else
 547     {
 548         md->fp_dhdl = fp_dhdl;
 549     }
 550     md->dhdl_derivatives = (ir->dhdl_derivatives==dhdlderivativesYES);
 551     return md;
 552 }
 553
 554 FILE *open_dhdl(const char *filename,const t_inputrec *ir,
 555                 const output_env_t oenv)
 556 {
 557     FILE *fp;
 558     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
 559     char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
 560     char **setname;
 561     char buf[STRLEN];
 562
 563     sprintf(label_x,"%s (%s)","Time",unit_time);
 564     if (ir->n_flambda == 0)
 565     {
 566         sprintf(title,"%s",dhdl);
 567         sprintf(label_y,"%s (%s %s)",
 568                 dhdl,unit_energy,"[\\lambda]\\S-1\\N");
 569     }
 570     else
 571     {
 572         sprintf(title,"%s, %s",dhdl,deltag);
 573         sprintf(label_y,"(%s)",unit_energy);
 574     }
 575     fp = gmx_fio_fopen(filename,"w+");
 576     xvgr_header(fp,title,label_x,label_y,exvggtXNY,oenv);
 577
 578     if (ir->delta_lambda == 0)
 579     {
 580         sprintf(buf,"T = %g (K), %s = %g",
 581                 ir->opts.ref_t[0],lambda,ir->init_lambda);
 582     }
 583     else
 584     {
 585         sprintf(buf,"T = %g (K)",
 586                 ir->opts.ref_t[0]);
 587     }
 588     xvgr_subtitle(fp,buf,oenv);
 589
 590     if (ir->n_flambda > 0)
 591     {
 592         int nsets,s,nsi=0;
 593         /* g_bar has to determine the lambda values used in this simulation
 594          * from this xvg legend.  */
 595         nsets = ( (ir->dhdl_derivatives==dhdlderivativesYES) ? 1 : 0) +
 596                   ir->n_flambda;
 597         snew(setname,nsets);
 598         if (ir->dhdl_derivatives == dhdlderivativesYES)
 599         {
 600             sprintf(buf,"%s %s %g",dhdl,lambda,ir->init_lambda);
 601             setname[nsi++] = strdup(buf);
 602         }
 603         for(s=0; s<ir->n_flambda; s++)
 604         {
 605             sprintf(buf,"%s %s %g",deltag,lambda,ir->flambda[s]);
 606             setname[nsi++] = strdup(buf);
 607         }
 608         xvgr_legend(fp,nsets,(const char**)setname,oenv);
 609
 610         for(s=0; s<nsets; s++)
 611         {
 612             sfree(setname[s]);
 613         }
 614         sfree(setname);
 615     }
 616
 617     return fp;
 618 }
 619
 620 static void copy_energy(t_mdebin *md, real e[],real ecpy[])
 621 {
 622     int i,j;
 623
 624     for(i=j=0; (i<F_NRE); i++)
 625         if (md->bEner[i])
 626             ecpy[j++] = e[i];
 627     if (j != md->f_nre)
 628         gmx_incons("Number of energy terms wrong");
 629 }
 630
 631 void upd_mdebin(t_mdebin *md, gmx_bool write_dhdl,
 632                 gmx_bool bSum,
 633                 double time,
 634                 real tmass,
 635                 gmx_enerdata_t *enerd,
 636                 t_state *state,
 637                 matrix  box,
 638                 tensor svir,
 639                 tensor fvir,
 640                 tensor vir,
 641                 tensor pres,
 642                 gmx_ekindata_t *ekind,
 643                 rvec mu_tot,
 644                 gmx_constr_t constr)
 645 {
 646     int    i,j,k,kk,m,n,gid;
 647     real   crmsd[2],tmp6[6];
 648     real   bs[NTRICLBOXS],vol,dens,pv,enthalpy;
 649     real   eee[egNR];
 650     real   ecopy[F_NRE];
 651     real   tmp;
 652     gmx_bool   bNoseHoover;
 653
 654     /* Do NOT use the box in the state variable, but the separate box provided
 655      * as an argument. This is because we sometimes need to write the box from
 656      * the last timestep to match the trajectory frames.
 657      */
 658     copy_energy(md, enerd->term,ecopy);
 659     add_ebin(md->ebin,md->ie,md->f_nre,ecopy,bSum);
 660     if (md->nCrmsd)
 661     {
 662         crmsd[0] = constr_rmsd(constr,FALSE);
 663         if (md->nCrmsd > 1)
 664         {
 665             crmsd[1] = constr_rmsd(constr,TRUE);
 666         }
 667         add_ebin(md->ebin,md->iconrmsd,md->nCrmsd,crmsd,FALSE);
 668     }
 669     if (md->bDynBox)
 670     {
 671         if(md->bTricl)
 672         {
 673             bs[0] = box[XX][XX];
 674             bs[1] = box[YY][XX];
 675             bs[2] = box[YY][YY];
 676             bs[3] = box[ZZ][XX];
 677             bs[4] = box[ZZ][YY];
 678             bs[5] = box[ZZ][ZZ];
 679         }
 680         else
 681         {
 682             bs[0] = box[XX][XX];
 683             bs[1] = box[YY][YY];
 684             bs[2] = box[ZZ][ZZ];
 685         }
 686         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 687         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 688
 689         /* This is pV (in kJ/mol).  The pressure is the reference pressure,
 690            not the instantaneous pressure */
 691         pv = 0;
 692         for (i=0;i<DIM;i++)
 693         {
 694             for (j=0;j<DIM;j++)
 695             {
 696                 if (i>j)
 697                 {
 698                     pv += box[i][j]*md->ref_p[i][j]/PRESFAC;
 699                 }
 700                 else
 701                 {
 702                     pv += box[j][i]*md->ref_p[j][i]/PRESFAC;
 703                 }
 704             }
 705         }
 706         add_ebin(md->ebin,md->ib   ,NBOXS,bs   ,bSum);
 707         add_ebin(md->ebin,md->ivol ,1    ,&vol ,bSum);
 708         add_ebin(md->ebin,md->idens,1    ,&dens,bSum);
 709         add_ebin(md->ebin,md->ipv  ,1    ,&pv  ,bSum);
 710         enthalpy = pv + enerd->term[F_ETOT];
 711         add_ebin(md->ebin,md->ienthalpy  ,1    ,&enthalpy  ,bSum);
 712     }
 713     if (md->bConstrVir)
 714     {
 715         add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
 716         add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
 717     }
 718     add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
 719     add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
 720     tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
 721     add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
 722     if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
 723     {
 724         tmp6[0] = state->boxv[XX][XX];
 725         tmp6[1] = state->boxv[YY][YY];
 726         tmp6[2] = state->boxv[ZZ][ZZ];
 727         tmp6[3] = state->boxv[YY][XX];
 728         tmp6[4] = state->boxv[ZZ][XX];
 729         tmp6[5] = state->boxv[ZZ][YY];
 730         add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
 731     }
 732     add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
 733     if (ekind && ekind->cosacc.cos_accel != 0)
 734     {
 735         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 736         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 737         add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
 738         /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
 739         tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
 740                  *vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
 741         add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
 742     }
 743     if (md->nE > 1)
 744     {
 745         n=0;
 746         for(i=0; (i<md->nEg); i++)
 747         {
 748             for(j=i; (j<md->nEg); j++)
 749             {
 750                 gid=GID(i,j,md->nEg);
 751                 for(k=kk=0; (k<egNR); k++)
 752                 {
 753                     if (md->bEInd[k])
 754                     {
 755                         eee[kk++] = enerd->grpp.ener[k][gid];
 756                     }
 757                 }
 758                 add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
 759                 n++;
 760             }
 761         }
 762     }
 763
 764     if (ekind)
 765     {
 766         for(i=0; (i<md->nTC); i++)
 767         {
 768             md->tmp_r[i] = ekind->tcstat[i].T;
 769         }
 770         add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
 771
 772         /* whether to print Nose-Hoover chains: */
 773         bNoseHoover = (getenv("GMX_NOSEHOOVER_CHAINS") != NULL);
 774
 775         if (md->etc == etcNOSEHOOVER)
 776         {
 777             if (bNoseHoover)
 778             {
 779                 if (md->bNHC_trotter)
 780                 {
 781                     for(i=0; (i<md->nTC); i++)
 782                     {
 783                         for (j=0;j<md->nNHC;j++)
 784                         {
 785                             k = i*md->nNHC+j;
 786                             md->tmp_r[2*k] = state->nosehoover_xi[k];
 787                             md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
 788                         }
 789                     }
 790                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
 791
 792                     if (md->bMTTK) {
 793                         for(i=0; (i<md->nTCP); i++)
 794                         {
 795                             for (j=0;j<md->nNHC;j++)
 796                             {
 797                                 k = i*md->nNHC+j;
 798                                 md->tmp_r[2*k] = state->nhpres_xi[k];
 799                                 md->tmp_r[2*k+1] = state->nhpres_vxi[k];
 800                             }
 801                         }
 802                         add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
 803                     }
 804
 805                 }
 806                 else
 807                 {
 808                     for(i=0; (i<md->nTC); i++)
 809                     {
 810                         md->tmp_r[2*i] = state->nosehoover_xi[i];
 811                         md->tmp_r[2*i+1] = state->nosehoover_vxi[i];
 812                     }
 813                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
 814                 }
 815             }
 816         }
 817         else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
 818                  md->etc == etcVRESCALE)
 819         {
 820             for(i=0; (i<md->nTC); i++)
 821             {
 822                 md->tmp_r[i] = ekind->tcstat[i].lambda;
 823             }
 824             add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
 825         }
 826     }
 827
 828     if (ekind && md->nU > 1)
 829     {
 830         for(i=0; (i<md->nU); i++)
 831         {
 832             copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
 833         }
 834         add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
 835     }
 836
 837     ebin_increase_count(md->ebin,bSum);
 838
 839     /* BAR + thermodynamic integration values */
 840     if (write_dhdl)
 841     {
 842         if (md->fp_dhdl)
 843         {
 844             fprintf(md->fp_dhdl,"%.4f", time);
 845
 846             if (md->dhdl_derivatives)
 847             {
 848                 fprintf(md->fp_dhdl," %g", enerd->term[F_DVDL]+
 849                                            enerd->term[F_DKDL]+
 850                                            enerd->term[F_DHDL_CON]);
 851             }
 852             for(i=1; i<enerd->n_lambda; i++)
 853             {
 854                 fprintf(md->fp_dhdl," %g",
 855                         enerd->enerpart_lambda[i]-enerd->enerpart_lambda[0]);
 856             }
 857             fprintf(md->fp_dhdl,"\n");
 858         }
 859         /* and the binary BAR output */
 860         if (md->dhc)
 861         {
 862             mde_delta_h_coll_add_dh(md->dhc,
 863                                     enerd->term[F_DVDL]+ enerd->term[F_DKDL]+
 864                                     enerd->term[F_DHDL_CON],
 865                                     enerd->enerpart_lambda, time,
 866                                     state->lambda);
 867         }
 868     }
 869 }
 870
 871 void upd_mdebin_step(t_mdebin *md)
 872 {
 873     ebin_increase_count(md->ebin,FALSE);
 874 }
 875
 876 static void npr(FILE *log,int n,char c)
 877 {
 878     for(; (n>0); n--) fprintf(log,"%c",c);
 879 }
 880
 881 static void pprint(FILE *log,const char *s,t_mdebin *md)
 882 {
 883     char CHAR='#';
 884     int  slen;
 885     char buf1[22],buf2[22];
 886
 887     slen = strlen(s);
 888     fprintf(log,"\t<======  ");
 889     npr(log,slen,CHAR);
 890     fprintf(log,"  ==>\n");
 891     fprintf(log,"\t<====  %s  ====>\n",s);
 892     fprintf(log,"\t<==  ");
 893     npr(log,slen,CHAR);
 894     fprintf(log,"  ======>\n\n");
 895
 896     fprintf(log,"\tStatistics over %s steps using %s frames\n",
 897             gmx_step_str(md->ebin->nsteps_sim,buf1),
 898             gmx_step_str(md->ebin->nsum_sim,buf2));
 899     fprintf(log,"\n");
 900 }
 901
 902 void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lamb)
 903 {
 904     char buf[22];
 905
 906     fprintf(log,"   %12s   %12s   %12s\n"
 907             "   %12s   %12.5f   %12.5f\n\n",
 908             "Step","Time","Lambda",gmx_step_str(steps,buf),time,lamb);
 909 }
 910
 911 void print_ebin(ener_file_t fp_ene,gmx_bool bEne,gmx_bool bDR,gmx_bool bOR,
 912                 FILE *log,
 913                 gmx_large_int_t step,double time,
 914                 int mode,gmx_bool bCompact,
 915                 t_mdebin *md,t_fcdata *fcd,
 916                 gmx_groups_t *groups,t_grpopts *opts)
 917 {
 918     /*static char **grpnms=NULL;*/
 919     char        buf[246];
 920     int         i,j,n,ni,nj,ndr,nor,b;
 921     int         ndisre=0;
 922     real        *disre_rm3tav, *disre_rt;
 923
 924     /* these are for the old-style blocks (1 subblock, only reals), because
 925        there can be only one per ID for these */
 926     int         nr[enxNR];
 927     int         id[enxNR];
 928     real        *block[enxNR];
 929
 930     /* temporary arrays for the lambda values to write out */
 931     double      enxlambda_data[2];
 932
 933     t_enxframe  fr;
 934
 935     switch (mode)
 936     {
 937         case eprNORMAL:
 938             init_enxframe(&fr);
 939             fr.t            = time;
 940             fr.step         = step;
 941             fr.nsteps       = md->ebin->nsteps;
 942             fr.dt           = md->delta_t;
 943             fr.nsum         = md->ebin->nsum;
 944             fr.nre          = (bEne) ? md->ebin->nener : 0;
 945             fr.ener         = md->ebin->e;
 946             ndisre          = bDR ? fcd->disres.npair : 0;
 947             disre_rm3tav    = fcd->disres.rm3tav;
 948             disre_rt        = fcd->disres.rt;
 949             /* Optional additional old-style (real-only) blocks. */
 950             for(i=0; i<enxNR; i++)
 951             {
 952                 nr[i] = 0;
 953             }
 954             if (fcd->orires.nr > 0 && bOR)
 955             {
 956                 diagonalize_orires_tensors(&(fcd->orires));
 957                 nr[enxOR]     = fcd->orires.nr;
 958                 block[enxOR]  = fcd->orires.otav;
 959                 id[enxOR]     = enxOR;
 960                 nr[enxORI]    = (fcd->orires.oinsl != fcd->orires.otav) ?
 961                           fcd->orires.nr : 0;
 962                 block[enxORI] = fcd->orires.oinsl;
 963                 id[enxORI]    = enxORI;
 964                 nr[enxORT]    = fcd->orires.nex*12;
 965                 block[enxORT] = fcd->orires.eig;
 966                 id[enxORT]    = enxORT;
 967             }
 968
 969             /* whether we are going to wrte anything out: */
 970             if (fr.nre || ndisre || nr[enxOR] || nr[enxORI])
 971             {
 972
 973                 /* the old-style blocks go first */
 974                 fr.nblock = 0;
 975                 for(i=0; i<enxNR; i++)
 976                 {
 977                     if (nr[i] > 0)
 978                     {
 979                         fr.nblock = i + 1;
 980                     }
 981                 }
 982                 add_blocks_enxframe(&fr, fr.nblock);
 983                 for(b=0;b<fr.nblock;b++)
 984                 {
 985                     add_subblocks_enxblock(&(fr.block[b]), 1);
 986                     fr.block[b].id=id[b];
 987                     fr.block[b].sub[0].nr = nr[b];
 988 #ifndef GMX_DOUBLE
 989                     fr.block[b].sub[0].type = xdr_datatype_float;
 990                     fr.block[b].sub[0].fval = block[b];
 991 #else
 992                     fr.block[b].sub[0].type = xdr_datatype_double;
 993                     fr.block[b].sub[0].dval = block[b];
 994 #endif
 995                 }
 996
 997                 /* check for disre block & fill it. */
 998                 if (ndisre>0)
 999                 {
1000                     int db = fr.nblock;
1001                     fr.nblock+=1;
1002                     add_blocks_enxframe(&fr, fr.nblock);
1003
1004                     add_subblocks_enxblock(&(fr.block[db]), 2);
1005                     fr.block[db].id=enxDISRE;
1006                     fr.block[db].sub[0].nr=ndisre;
1007                     fr.block[db].sub[1].nr=ndisre;
1008 #ifndef GMX_DOUBLE
1009                     fr.block[db].sub[0].type=xdr_datatype_float;
1010                     fr.block[db].sub[1].type=xdr_datatype_float;
1011                     fr.block[db].sub[0].fval=disre_rt;
1012                     fr.block[db].sub[1].fval=disre_rm3tav;
1013 #else
1014                     fr.block[db].sub[0].type=xdr_datatype_double;
1015                     fr.block[db].sub[1].type=xdr_datatype_double;
1016                     fr.block[db].sub[0].dval=disre_rt;
1017                     fr.block[db].sub[1].dval=disre_rm3tav;
1018 #endif
1019                 }
1020                 /* here we can put new-style blocks */
1021
1022                 /* Free energy perturbation blocks */
1023                 if (md->dhc)
1024                 {
1025                     mde_delta_h_coll_handle_block(md->dhc, &fr, fr.nblock);
1026                 }
1027
1028                 /* do the actual I/O */
1029                 do_enx(fp_ene,&fr);
1030                 gmx_fio_check_file_position(enx_file_pointer(fp_ene));
1031                 if (fr.nre)
1032                 {
1033                     /* We have stored the sums, so reset the sum history */
1034                     reset_ebin_sums(md->ebin);
1035                 }
1036
1037                 /* we can now free & reset the data in the blocks */
1038                 if (md->dhc)
1039                     mde_delta_h_coll_reset(md->dhc);
1040             }
1041             free_enxframe(&fr);
1042             break;
1043         case eprAVER:
1044             if (log)
1045             {
1046                 pprint(log,"A V E R A G E S",md);
1047             }
1048             break;
1049         case eprRMS:
1050             if (log)
1051             {
1052                 pprint(log,"R M S - F L U C T U A T I O N S",md);
1053             }
1054             break;
1055         default:
1056             gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
1057     }
1058
1059     if (log)
1060     {
1061         for(i=0;i<opts->ngtc;i++)
1062         {
1063             if(opts->annealing[i]!=eannNO)
1064             {
1065                 fprintf(log,"Current ref_t for group %s: %8.1f\n",
1066                         *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
1067                         opts->ref_t[i]);
1068             }
1069         }
1070         if (mode==eprNORMAL && fcd->orires.nr>0)
1071         {
1072             print_orires_log(log,&(fcd->orires));
1073         }
1074         fprintf(log,"   Energies (%s)\n",unit_energy);
1075         pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
1076         fprintf(log,"\n");
1077
1078         if (!bCompact)
1079         {
1080             if (md->bDynBox)
1081             {
1082                 pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
1083                         mode,TRUE);
1084                 fprintf(log,"\n");
1085             }
1086             if (md->bConstrVir)
1087             {
1088                 fprintf(log,"   Constraint Virial (%s)\n",unit_energy);
1089                 pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
1090                 fprintf(log,"\n");
1091                 fprintf(log,"   Force Virial (%s)\n",unit_energy);
1092                 pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
1093                 fprintf(log,"\n");
1094             }
1095             fprintf(log,"   Total Virial (%s)\n",unit_energy);
1096             pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
1097             fprintf(log,"\n");
1098             fprintf(log,"   Pressure (%s)\n",unit_pres_bar);
1099             pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
1100             fprintf(log,"\n");
1101             fprintf(log,"   Total Dipole (%s)\n",unit_dipole_D);
1102             pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
1103             fprintf(log,"\n");
1104
1105             if (md->nE > 1)
1106             {
1107                 if (md->print_grpnms==NULL)
1108                 {
1109                     snew(md->print_grpnms,md->nE);
1110                     n=0;
1111                     for(i=0; (i<md->nEg); i++)
1112                     {
1113                         ni=groups->grps[egcENER].nm_ind[i];
1114                         for(j=i; (j<md->nEg); j++)
1115                         {
1116                             nj=groups->grps[egcENER].nm_ind[j];
1117                             sprintf(buf,"%s-%s",*(groups->grpname[ni]),
1118                                     *(groups->grpname[nj]));
1119                             md->print_grpnms[n++]=strdup(buf);
1120                         }
1121                     }
1122                 }
1123                 sprintf(buf,"Epot (%s)",unit_energy);
1124                 fprintf(log,"%15s   ",buf);
1125                 for(i=0; (i<egNR); i++)
1126                 {
1127                     if (md->bEInd[i])
1128                     {
1129                         fprintf(log,"%12s   ",egrp_nm[i]);
1130                     }
1131                 }
1132                 fprintf(log,"\n");
1133                 for(i=0; (i<md->nE); i++)
1134                 {
1135                     fprintf(log,"%15s",md->print_grpnms[i]);
1136                     pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
1137                             FALSE);
1138                 }
1139                 fprintf(log,"\n");
1140             }
1141             if (md->nTC > 1)
1142             {
1143                 pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
1144                 fprintf(log,"\n");
1145             }
1146             if (md->nU > 1)
1147             {
1148                 fprintf(log,"%15s   %12s   %12s   %12s\n",
1149                         "Group","Ux","Uy","Uz");
1150                 for(i=0; (i<md->nU); i++)
1151                 {
1152                     ni=groups->grps[egcACC].nm_ind[i];
1153                     fprintf(log,"%15s",*groups->grpname[ni]);
1154                     pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
1155                 }
1156                 fprintf(log,"\n");
1157             }
1158         }
1159     }
1160
1161 }
1162
1163 void update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
1164 {
1165     int i;
1166
1167     enerhist->nsteps     = mdebin->ebin->nsteps;
1168     enerhist->nsum       = mdebin->ebin->nsum;
1169     enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
1170     enerhist->nsum_sim   = mdebin->ebin->nsum_sim;
1171     enerhist->nener      = mdebin->ebin->nener;
1172
1173     if (mdebin->ebin->nsum > 0)
1174     {
1175         /* Check if we need to allocate first */
1176         if(enerhist->ener_ave == NULL)
1177         {
1178             snew(enerhist->ener_ave,enerhist->nener);
1179             snew(enerhist->ener_sum,enerhist->nener);
1180         }
1181
1182         for(i=0;i<enerhist->nener;i++)
1183         {
1184             enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
1185             enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
1186         }
1187     }
1188
1189     if (mdebin->ebin->nsum_sim > 0)
1190     {
1191         /* Check if we need to allocate first */
1192         if(enerhist->ener_sum_sim == NULL)
1193         {
1194             snew(enerhist->ener_sum_sim,enerhist->nener);
1195         }
1196
1197         for(i=0;i<enerhist->nener;i++)
1198         {
1199             enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
1200         }
1201     }
1202     if (mdebin->dhc)
1203     {
1204         mde_delta_h_coll_update_energyhistory(mdebin->dhc, enerhist);
1205     }
1206 }
1207
1208 void restore_energyhistory_from_state(t_mdebin * mdebin,
1209                                       energyhistory_t * enerhist)
1210 {
1211     int i;
1212
1213     if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
1214         mdebin->ebin->nener != enerhist->nener)
1215     {
1216         gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1217                   mdebin->ebin->nener,enerhist->nener);
1218     }
1219
1220     mdebin->ebin->nsteps     = enerhist->nsteps;
1221     mdebin->ebin->nsum       = enerhist->nsum;
1222     mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
1223     mdebin->ebin->nsum_sim   = enerhist->nsum_sim;
1224
1225     for(i=0; i<mdebin->ebin->nener; i++)
1226     {
1227         mdebin->ebin->e[i].eav  =
1228                   (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
1229         mdebin->ebin->e[i].esum =
1230                   (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
1231         mdebin->ebin->e_sim[i].esum =
1232                   (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
1233     }
1234     if (mdebin->dhc)
1235     {
1236         mde_delta_h_coll_restore_energyhistory(mdebin->dhc, enerhist);
1237     }
1238 }