src/gromacs/mdlib/mdebin.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <string.h>
  41 #include <float.h>
  42 #include "typedefs.h"
  43 #include "string2.h"
  44 #include "mdebin.h"
  45 #include "smalloc.h"
  46 #include "physics.h"
  47 #include "enxio.h"
  48 #include "vec.h"
  49 #include "disre.h"
  50 #include "main.h"
  51 #include "network.h"
  52 #include "names.h"
  53 #include "orires.h"
  54 #include "constr.h"
  55 #include "mtop_util.h"
  56 #include "xvgr.h"
  57 #include "gmxfio.h"
  58 #include "macros.h"
  59 #include "mdrun.h"
  60 #include "mdebin_bar.h"
  61
  62
  63 static const char *conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
  64
  65 static const char *boxs_nm[] = { "Box-X", "Box-Y", "Box-Z" };
  66
  67 static const char *tricl_boxs_nm[] = {
  68     "Box-XX", "Box-YY", "Box-ZZ",
  69     "Box-YX", "Box-ZX", "Box-ZY"
  70 };
  71
  72 static const char *vol_nm[] = { "Volume" };
  73
  74 static const char *dens_nm[] = {"Density" };
  75
  76 static const char *pv_nm[] = {"pV" };
  77
  78 static const char *enthalpy_nm[] = {"Enthalpy" };
  79
  80 static const char *boxvel_nm[] = {
  81     "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
  82     "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
  83 };
  84
  85 #define NBOXS asize(boxs_nm)
  86 #define NTRICLBOXS asize(tricl_boxs_nm)
  87
  88 t_mdebin *init_mdebin(ener_file_t fp_ene,
  89                       const gmx_mtop_t *mtop,
  90                       const t_inputrec *ir,
  91                       FILE *fp_dhdl)
  92 {
  93     const char *ener_nm[F_NRE];
  94     static const char *vir_nm[] = {
  95         "Vir-XX", "Vir-XY", "Vir-XZ",
  96         "Vir-YX", "Vir-YY", "Vir-YZ",
  97         "Vir-ZX", "Vir-ZY", "Vir-ZZ"
  98     };
  99     static const char *sv_nm[] = {
 100         "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
 101         "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
 102         "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
 103     };
 104     static const char *fv_nm[] = {
 105         "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
 106         "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
 107         "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
 108     };
 109     static const char *pres_nm[] = {
 110         "Pres-XX","Pres-XY","Pres-XZ",
 111         "Pres-YX","Pres-YY","Pres-YZ",
 112         "Pres-ZX","Pres-ZY","Pres-ZZ"
 113     };
 114     static const char *surft_nm[] = {
 115         "#Surf*SurfTen"
 116     };
 117     static const char *mu_nm[] = {
 118         "Mu-X", "Mu-Y", "Mu-Z"
 119     };
 120     static const char *vcos_nm[] = {
 121         "2CosZ*Vel-X"
 122     };
 123     static const char *visc_nm[] = {
 124         "1/Viscosity"
 125     };
 126     static const char *baro_nm[] = {
 127         "Barostat"
 128     };
 129
 130     char     **grpnms;
 131     const gmx_groups_t *groups;
 132     char     **gnm;
 133     char     buf[256];
 134     const char     *bufi;
 135     t_mdebin *md;
 136     int      i,j,ni,nj,n,nh,k,kk,ncon,nset;
 137     gmx_bool     bBHAM,bNoseHoover,b14;
 138
 139     snew(md,1);
 140
 141     md->bVir=TRUE;
 142     md->bPress=TRUE;
 143     md->bSurft=TRUE;
 144     md->bMu=TRUE;
 145
 146     if (EI_DYNAMICS(ir->eI))
 147     {
 148         md->delta_t = ir->delta_t;
 149     }
 150     else
 151     {
 152         md->delta_t = 0;
 153     }
 154
 155     groups = &mtop->groups;
 156
 157     bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
 158     b14   = (gmx_mtop_ftype_count(mtop,F_LJ14) > 0 ||
 159              gmx_mtop_ftype_count(mtop,F_LJC14_Q) > 0);
 160
 161     ncon = gmx_mtop_ftype_count(mtop,F_CONSTR);
 162     nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
 163     md->bConstr    = (ncon > 0 || nset > 0);
 164     md->bConstrVir = FALSE;
 165     if (md->bConstr) {
 166         if (ncon > 0 && ir->eConstrAlg == econtLINCS) {
 167             if (ir->eI == eiSD2)
 168                 md->nCrmsd = 2;
 169             else
 170                 md->nCrmsd = 1;
 171         }
 172         md->bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL);
 173     } else {
 174         md->nCrmsd = 0;
 175     }
 176
 177     /* Energy monitoring */
 178     for(i=0;i<egNR;i++)
 179     {
 180         md->bEInd[i]=FALSE;
 181     }
 182
 183     /* Even though the OpenMM build has moved to contrib, it's not
 184      * practical to move/remove this code fragment, because of the
 185      * fundamental mess that is the GROMACS library structure. */
 186 #ifndef GMX_OPENMM
 187     for(i=0; i<F_NRE; i++)
 188     {
 189         md->bEner[i] = FALSE;
 190         if (i == F_LJ)
 191             md->bEner[i] = !bBHAM;
 192         else if (i == F_BHAM)
 193             md->bEner[i] = bBHAM;
 194         else if (i == F_EQM)
 195             md->bEner[i] = ir->bQMMM;
 196         else if (i == F_COUL_LR)
 197             md->bEner[i] = (ir->rcoulomb > ir->rlist);
 198         else if (i == F_LJ_LR)
 199             md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist);
 200         else if (i == F_BHAM_LR)
 201             md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist);
 202         else if (i == F_RF_EXCL)
 203             md->bEner[i] = (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC && ir->cutoff_scheme == ecutsGROUP);
 204         else if (i == F_COUL_RECIP)
 205             md->bEner[i] = EEL_FULL(ir->coulombtype);
 206         else if (i == F_LJ14)
 207             md->bEner[i] = b14;
 208         else if (i == F_COUL14)
 209             md->bEner[i] = b14;
 210         else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
 211             md->bEner[i] = FALSE;
 212         else if ((i == F_DVDL_COUL && ir->fepvals->separate_dvdl[efptCOUL]) ||
 213                  (i == F_DVDL_VDW  && ir->fepvals->separate_dvdl[efptVDW]) ||
 214                  (i == F_DVDL_BONDED && ir->fepvals->separate_dvdl[efptBONDED]) ||
 215                  (i == F_DVDL_RESTRAINT && ir->fepvals->separate_dvdl[efptRESTRAINT]) ||
 216                  (i == F_DKDL && ir->fepvals->separate_dvdl[efptMASS]) ||
 217                  (i == F_DVDL && ir->fepvals->separate_dvdl[efptFEP]))
 218             md->bEner[i] = (ir->efep != efepNO);
 219         else if ((interaction_function[i].flags & IF_VSITE) ||
 220                  (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
 221             md->bEner[i] = FALSE;
 222         else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES)  || (i==F_EQM))
 223             md->bEner[i] = TRUE;
 224         else if ((i == F_GBPOL) && ir->implicit_solvent==eisGBSA)
 225             md->bEner[i] = TRUE;
 226         else if ((i == F_NPSOLVATION) && ir->implicit_solvent==eisGBSA && (ir->sa_algorithm != esaNO))
 227             md->bEner[i] = TRUE;
 228         else if ((i == F_GB12) || (i == F_GB13) || (i == F_GB14))
 229             md->bEner[i] = FALSE;
 230         else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
 231             md->bEner[i] = EI_DYNAMICS(ir->eI);
 232         else if (i == F_DISPCORR || i == F_PDISPCORR)
 233             md->bEner[i] = (ir->eDispCorr != edispcNO);
 234         else if (i == F_DISRESVIOL)
 235             md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_DISRES) > 0);
 236         else if (i == F_ORIRESDEV)
 237             md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0);
 238         else if (i == F_CONNBONDS)
 239             md->bEner[i] = FALSE;
 240         else if (i == F_COM_PULL)
 241             md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F || ir->bRot);
 242         else if (i == F_ECONSERVED)
 243             md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) &&
 244                             (ir->epc == epcNO || ir->epc==epcMTTK));
 245         else
 246             md->bEner[i] = (gmx_mtop_ftype_count(mtop,i) > 0);
 247     }
 248 #else
 249     /* OpenMM always produces only the following 4 energy terms */
 250     md->bEner[F_EPOT] = TRUE;
 251     md->bEner[F_EKIN] = TRUE;
 252     md->bEner[F_ETOT] = TRUE;
 253     md->bEner[F_TEMP] = TRUE;
 254 #endif
 255
 256     /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
 257     if (ir->bAdress && !debug) {
 258         for (i = 0; i < F_NRE; i++) {
 259             md->bEner[i] = FALSE;
 260             if(i == F_EKIN){ md->bEner[i] = TRUE;}
 261             if(i == F_TEMP){ md->bEner[i] = TRUE;}
 262         }
 263         md->bVir=FALSE;
 264         md->bPress=FALSE;
 265         md->bSurft=FALSE;
 266         md->bMu=FALSE;
 267     }
 268
 269     md->f_nre=0;
 270     for(i=0; i<F_NRE; i++)
 271     {
 272         if (md->bEner[i])
 273         {
 274             ener_nm[md->f_nre]=interaction_function[i].longname;
 275             md->f_nre++;
 276         }
 277     }
 278
 279     md->epc = ir->epc;
 280     md->bDiagPres = !TRICLINIC(ir->ref_p);
 281     md->ref_p = (ir->ref_p[XX][XX]+ir->ref_p[YY][YY]+ir->ref_p[ZZ][ZZ])/DIM;
 282     md->bTricl = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
 283     md->bDynBox = DYNAMIC_BOX(*ir);
 284     md->etc = ir->etc;
 285     md->bNHC_trotter = IR_NVT_TROTTER(ir);
 286     md->bPrintNHChains = ir-> bPrintNHChains;
 287     md->bMTTK = (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir));
 288     md->bMu = NEED_MUTOT(*ir);
 289
 290     md->ebin  = mk_ebin();
 291     /* Pass NULL for unit to let get_ebin_space determine the units
 292      * for interaction_function[i].longname
 293      */
 294     md->ie    = get_ebin_space(md->ebin,md->f_nre,ener_nm,NULL);
 295     if (md->nCrmsd)
 296     {
 297         /* This should be called directly after the call for md->ie,
 298          * such that md->iconrmsd follows directly in the list.
 299          */
 300         md->iconrmsd = get_ebin_space(md->ebin,md->nCrmsd,conrmsd_nm,"");
 301     }
 302     if (md->bDynBox)
 303     {
 304         md->ib    = get_ebin_space(md->ebin,
 305                                    md->bTricl ? NTRICLBOXS : NBOXS,
 306                                    md->bTricl ? tricl_boxs_nm : boxs_nm,
 307                                    unit_length);
 308         md->ivol  = get_ebin_space(md->ebin, 1, vol_nm,  unit_volume);
 309         md->idens = get_ebin_space(md->ebin, 1, dens_nm, unit_density_SI);
 310         if (md->bDiagPres)
 311         {
 312             md->ipv   = get_ebin_space(md->ebin, 1, pv_nm,   unit_energy);
 313             md->ienthalpy = get_ebin_space(md->ebin, 1, enthalpy_nm,   unit_energy);
 314         }
 315     }
 316     if (md->bConstrVir)
 317     {
 318         md->isvir = get_ebin_space(md->ebin,asize(sv_nm),sv_nm,unit_energy);
 319         md->ifvir = get_ebin_space(md->ebin,asize(fv_nm),fv_nm,unit_energy);
 320     }
 321     if (md->bVir)
 322         md->ivir   = get_ebin_space(md->ebin,asize(vir_nm),vir_nm,unit_energy);
 323     if (md->bPress)
 324         md->ipres  = get_ebin_space(md->ebin,asize(pres_nm),pres_nm,unit_pres_bar);
 325     if (md->bSurft)
 326         md->isurft = get_ebin_space(md->ebin,asize(surft_nm),surft_nm,
 327                                 unit_surft_bar);
 328     if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
 329     {
 330         md->ipc = get_ebin_space(md->ebin,md->bTricl ? 6 : 3,
 331                                  boxvel_nm,unit_vel);
 332     }
 333     if (md->bMu)
 334     {
 335         md->imu    = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
 336     }
 337     if (ir->cos_accel != 0)
 338     {
 339         md->ivcos = get_ebin_space(md->ebin,asize(vcos_nm),vcos_nm,unit_vel);
 340         md->ivisc = get_ebin_space(md->ebin,asize(visc_nm),visc_nm,
 341                                    unit_invvisc_SI);
 342     }
 343
 344     /* Energy monitoring */
 345     for(i=0;i<egNR;i++)
 346     {
 347         md->bEInd[i] = FALSE;
 348     }
 349     md->bEInd[egCOULSR] = TRUE;
 350     md->bEInd[egLJSR  ] = TRUE;
 351
 352     if (ir->rcoulomb > ir->rlist)
 353     {
 354         md->bEInd[egCOULLR] = TRUE;
 355     }
 356     if (!bBHAM)
 357     {
 358         if (ir->rvdw > ir->rlist)
 359         {
 360             md->bEInd[egLJLR]   = TRUE;
 361         }
 362     }
 363     else
 364     {
 365         md->bEInd[egLJSR]   = FALSE;
 366         md->bEInd[egBHAMSR] = TRUE;
 367         if (ir->rvdw > ir->rlist)
 368         {
 369             md->bEInd[egBHAMLR]   = TRUE;
 370         }
 371     }
 372     if (b14)
 373     {
 374         md->bEInd[egLJ14] = TRUE;
 375         md->bEInd[egCOUL14] = TRUE;
 376     }
 377     md->nEc=0;
 378     for(i=0; (i<egNR); i++)
 379     {
 380         if (md->bEInd[i])
 381         {
 382             md->nEc++;
 383         }
 384     }
 385
 386     n=groups->grps[egcENER].nr;
 387     /* for adress simulations, most energy terms are not meaningfull, and thus disabled*/
 388     if (!ir->bAdress){
 389         /*standard simulation*/
 390         md->nEg=n;
 391         md->nE=(n*(n+1))/2;
 392     }
 393     else if (!debug) {
 394         /*AdResS simulation*/
 395        md->nU=0;
 396        md->nEg=0;
 397        md->nE=0;
 398        md->nEc=0;
 399        md->isvir=FALSE;
 400     }
 401     snew(md->igrp,md->nE);
 402     if (md->nE > 1)
 403     {
 404         n=0;
 405         snew(gnm,md->nEc);
 406         for(k=0; (k<md->nEc); k++)
 407         {
 408             snew(gnm[k],STRLEN);
 409         }
 410         for(i=0; (i<groups->grps[egcENER].nr); i++)
 411         {
 412             ni=groups->grps[egcENER].nm_ind[i];
 413             for(j=i; (j<groups->grps[egcENER].nr); j++)
 414             {
 415                 nj=groups->grps[egcENER].nm_ind[j];
 416                 for(k=kk=0; (k<egNR); k++)
 417                 {
 418                     if (md->bEInd[k])
 419                     {
 420                         sprintf(gnm[kk],"%s:%s-%s",egrp_nm[k],
 421                                 *(groups->grpname[ni]),*(groups->grpname[nj]));
 422                         kk++;
 423                     }
 424                 }
 425                 md->igrp[n]=get_ebin_space(md->ebin,md->nEc,
 426                                            (const char **)gnm,unit_energy);
 427                 n++;
 428             }
 429         }
 430         for(k=0; (k<md->nEc); k++)
 431         {
 432             sfree(gnm[k]);
 433         }
 434         sfree(gnm);
 435
 436         if (n != md->nE)
 437         {
 438             gmx_incons("Number of energy terms wrong");
 439         }
 440     }
 441
 442     md->nTC=groups->grps[egcTC].nr;
 443     md->nNHC = ir->opts.nhchainlength; /* shorthand for number of NH chains */
 444     if (md->bMTTK)
 445     {
 446         md->nTCP = 1;  /* assume only one possible coupling system for barostat
 447                           for now */
 448     }
 449     else
 450     {
 451         md->nTCP = 0;
 452     }
 453     if (md->etc == etcNOSEHOOVER)
 454     {
 455         if (md->bNHC_trotter)
 456         {
 457             md->mde_n = 2*md->nNHC*md->nTC;
 458         }
 459         else
 460         {
 461             md->mde_n = 2*md->nTC;
 462         }
 463         if (md->epc == epcMTTK)
 464         {
 465             md->mdeb_n = 2*md->nNHC*md->nTCP;
 466         }
 467     } else {
 468         md->mde_n = md->nTC;
 469         md->mdeb_n = 0;
 470     }
 471
 472     snew(md->tmp_r,md->mde_n);
 473     snew(md->tmp_v,md->mde_n);
 474     snew(md->grpnms,md->mde_n);
 475     grpnms = md->grpnms;
 476
 477     for(i=0; (i<md->nTC); i++)
 478     {
 479         ni=groups->grps[egcTC].nm_ind[i];
 480         sprintf(buf,"T-%s",*(groups->grpname[ni]));
 481         grpnms[i]=strdup(buf);
 482     }
 483     md->itemp=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,
 484                              unit_temp_K);
 485
 486     if (md->etc == etcNOSEHOOVER)
 487     {
 488         if (md->bPrintNHChains)
 489         {
 490             if (md->bNHC_trotter)
 491             {
 492                 for(i=0; (i<md->nTC); i++)
 493                 {
 494                     ni=groups->grps[egcTC].nm_ind[i];
 495                     bufi = *(groups->grpname[ni]);
 496                     for(j=0; (j<md->nNHC); j++)
 497                     {
 498                         sprintf(buf,"Xi-%d-%s",j,bufi);
 499                         grpnms[2*(i*md->nNHC+j)]=strdup(buf);
 500                         sprintf(buf,"vXi-%d-%s",j,bufi);
 501                         grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
 502                     }
 503                 }
 504                 md->itc=get_ebin_space(md->ebin,md->mde_n,
 505                                        (const char **)grpnms,unit_invtime);
 506                 if (md->bMTTK)
 507                 {
 508                     for(i=0; (i<md->nTCP); i++)
 509                     {
 510                         bufi = baro_nm[0];  /* All barostat DOF's together for now. */
 511                         for(j=0; (j<md->nNHC); j++)
 512                         {
 513                             sprintf(buf,"Xi-%d-%s",j,bufi);
 514                             grpnms[2*(i*md->nNHC+j)]=strdup(buf);
 515                             sprintf(buf,"vXi-%d-%s",j,bufi);
 516                             grpnms[2*(i*md->nNHC+j)+1]=strdup(buf);
 517                         }
 518                     }
 519                     md->itcb=get_ebin_space(md->ebin,md->mdeb_n,
 520                                             (const char **)grpnms,unit_invtime);
 521                 }
 522             }
 523             else
 524             {
 525                 for(i=0; (i<md->nTC); i++)
 526                 {
 527                     ni=groups->grps[egcTC].nm_ind[i];
 528                     bufi = *(groups->grpname[ni]);
 529                     sprintf(buf,"Xi-%s",bufi);
 530                     grpnms[2*i]=strdup(buf);
 531                     sprintf(buf,"vXi-%s",bufi);
 532                     grpnms[2*i+1]=strdup(buf);
 533                 }
 534                 md->itc=get_ebin_space(md->ebin,md->mde_n,
 535                                        (const char **)grpnms,unit_invtime);
 536             }
 537         }
 538     }
 539     else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
 540              md->etc == etcVRESCALE)
 541     {
 542         for(i=0; (i<md->nTC); i++)
 543         {
 544             ni=groups->grps[egcTC].nm_ind[i];
 545             sprintf(buf,"Lamb-%s",*(groups->grpname[ni]));
 546             grpnms[i]=strdup(buf);
 547         }
 548         md->itc=get_ebin_space(md->ebin,md->mde_n,(const char **)grpnms,"");
 549     }
 550
 551     sfree(grpnms);
 552
 553
 554     md->nU=groups->grps[egcACC].nr;
 555     if (md->nU > 1)
 556     {
 557         snew(grpnms,3*md->nU);
 558         for(i=0; (i<md->nU); i++)
 559         {
 560             ni=groups->grps[egcACC].nm_ind[i];
 561             sprintf(buf,"Ux-%s",*(groups->grpname[ni]));
 562             grpnms[3*i+XX]=strdup(buf);
 563             sprintf(buf,"Uy-%s",*(groups->grpname[ni]));
 564             grpnms[3*i+YY]=strdup(buf);
 565             sprintf(buf,"Uz-%s",*(groups->grpname[ni]));
 566             grpnms[3*i+ZZ]=strdup(buf);
 567         }
 568         md->iu=get_ebin_space(md->ebin,3*md->nU,(const char **)grpnms,unit_vel);
 569         sfree(grpnms);
 570     }
 571
 572     if ( fp_ene )
 573     {
 574         do_enxnms(fp_ene,&md->ebin->nener,&md->ebin->enm);
 575     }
 576
 577     md->print_grpnms=NULL;
 578
 579     /* check whether we're going to write dh histograms */
 580     md->dhc=NULL;
 581     if (ir->fepvals->separate_dhdl_file == esepdhdlfileNO )
 582     {
 583         /* Currently dh histograms are only written with dynamics */
 584         if (EI_DYNAMICS(ir->eI))
 585         {
 586             snew(md->dhc, 1);
 587
 588             mde_delta_h_coll_init(md->dhc, ir);
 589         }
 590         md->fp_dhdl = NULL;
 591         snew(md->dE,ir->fepvals->n_lambda);
 592     }
 593     else
 594     {
 595         md->fp_dhdl = fp_dhdl;
 596         snew(md->dE,ir->fepvals->n_lambda);
 597     }
 598     if (ir->bSimTemp) {
 599         int i;
 600         snew(md->temperatures,ir->fepvals->n_lambda);
 601         for (i=0;i<ir->fepvals->n_lambda;i++)
 602         {
 603             md->temperatures[i] = ir->simtempvals->temperatures[i];
 604         }
 605     }
 606     return md;
 607 }
 608
 609 /* print a lambda vector to a string
 610    fep = the inputrec's FEP input data
 611    i = the index of the lambda vector
 612    get_native_lambda = whether to print the native lambda
 613    get_names = whether to print the names rather than the values
 614    str = the pre-allocated string buffer to print to. */
 615 static void print_lambda_vector(t_lambda *fep, int i,
 616                                 gmx_bool get_native_lambda, gmx_bool get_names,
 617                                 char *str)
 618 {
 619     size_t nps=0, np;
 620     int j,k=0;
 621     int Nsep=0;
 622
 623     for (j=0;j<efptNR;j++)
 624     {
 625         if (fep->separate_dvdl[j])
 626             Nsep ++;
 627     }
 628     str[0]=0; /* reset the string */
 629     if (Nsep > 1)
 630     {
 631         str += sprintf(str, "("); /* set the opening parenthesis*/
 632     }
 633     for (j=0;j<efptNR;j++)
 634     {
 635         if (fep->separate_dvdl[j])
 636         {
 637             double lam;
 638             if (!get_names)
 639             {
 640                 if (get_native_lambda && fep->init_lambda >= 0)
 641                 {
 642                     str += sprintf(str,"%.4f", fep->init_lambda);
 643                 }
 644                 else
 645                 {
 646                     str += sprintf(str,"%.4f", fep->all_lambda[j][i]);
 647                 }
 648             }
 649             else
 650             {
 651                 str += sprintf(str,"%s", efpt_singular_names[j]);
 652             }
 653             /* print comma for the next item */
 654             if (k<Nsep-1)
 655             {
 656                 str += sprintf(str,", ");
 657             }
 658             k++;
 659         }
 660     }
 661     if (Nsep > 1)
 662     {
 663         /* and add the closing parenthesis */
 664         str += sprintf(str, ")");
 665     }
 666 }
 667
 668
 669 extern FILE *open_dhdl(const char *filename,const t_inputrec *ir,
 670                        const output_env_t oenv)
 671 {
 672     FILE *fp;
 673     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda",
 674         *lambdastate="\\lambda state",*remain="remaining";
 675     char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
 676     int  i,np,nps,nsets,nsets_de,nsetsbegin;
 677     int n_lambda_terms=0;
 678     t_lambda *fep=ir->fepvals; /* for simplicity */
 679     t_expanded *expand=ir->expandedvals;
 680     char **setname;
 681     char buf[STRLEN], lambda_vec_str[STRLEN], lambda_name_str[STRLEN];
 682     int bufplace=0;
 683
 684     int nsets_dhdl = 0;
 685     int s = 0;
 686     int nsetsextend;
 687     gmx_bool write_pV = FALSE;
 688
 689     /* count the number of different lambda terms */
 690     for (i=0;i<efptNR;i++)
 691     {
 692         if (fep->separate_dvdl[i])
 693         {
 694             n_lambda_terms++;
 695         }
 696     }
 697
 698     if (fep->n_lambda == 0)
 699     {
 700         sprintf(title,"%s",dhdl);
 701         sprintf(label_x,"Time (ps)");
 702         sprintf(label_y,"%s (%s %s)",
 703                 dhdl,unit_energy,"[\\lambda]\\S-1\\N");
 704     }
 705     else
 706     {
 707         sprintf(title,"%s and %s",dhdl,deltag);
 708         sprintf(label_x,"Time (ps)");
 709         sprintf(label_y,"%s and %s (%s %s)",
 710                 dhdl,deltag,unit_energy,"[\\8l\\4]\\S-1\\N");
 711     }
 712     fp = gmx_fio_fopen(filename,"w+");
 713     xvgr_header(fp,title,label_x,label_y,exvggtXNY,oenv);
 714
 715     if (!(ir->bSimTemp))
 716     {
 717         bufplace = sprintf(buf,"T = %g (K) ",
 718                 ir->opts.ref_t[0]);
 719     }
 720     if (ir->efep != efepSLOWGROWTH)
 721     {
 722         if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))
 723         {
 724             /* compatibility output */
 725             sprintf(&(buf[bufplace]),"%s = %.4f", lambda,fep->init_lambda);
 726         }
 727         else
 728         {
 729             print_lambda_vector(fep, fep->init_fep_state, TRUE, FALSE,
 730                                 lambda_vec_str);
 731             print_lambda_vector(fep, fep->init_fep_state, TRUE, TRUE,
 732                                 lambda_name_str);
 733             sprintf(&(buf[bufplace]),"%s %d: %s = %s",
 734                     lambdastate,fep->init_fep_state,
 735                     lambda_name_str, lambda_vec_str);
 736         }
 737     }
 738     xvgr_subtitle(fp,buf,oenv);
 739
 740
 741     nsets_dhdl=0;
 742     if (fep->dhdl_derivatives == edhdlderivativesYES)
 743     {
 744         nsets_dhdl = n_lambda_terms;
 745     }
 746     /* count the number of delta_g states */
 747     nsets_de = fep->lambda_stop_n - fep->lambda_start_n;
 748
 749     nsets = nsets_dhdl + nsets_de; /* dhdl + fep differences */
 750
 751     if (fep->n_lambda>0 && (expand->elmcmove > elmcmoveNO))
 752     {
 753         nsets += 1;   /*add fep state for expanded ensemble */
 754     }
 755
 756     if (fep->bPrintEnergy)
 757     {
 758         nsets += 1;  /* add energy to the dhdl as well */
 759     }
 760
 761     nsetsextend = nsets;
 762     if ((ir->epc!=epcNO) && (fep->n_lambda>0) && (fep->init_lambda < 0))
 763     {
 764         nsetsextend += 1; /* for PV term, other terms possible if required for
 765                              the reduced potential (only needed with foreign
 766                              lambda, and only output when init_lambda is not
 767                              set in order to maintain compatibility of the
 768                              dhdl.xvg file) */
 769         write_pV = TRUE;
 770     }
 771     snew(setname,nsetsextend);
 772
 773     if (expand->elmcmove > elmcmoveNO)
 774     {
 775         /* state for the fep_vals, if we have alchemical sampling */
 776         sprintf(buf,"%s","Thermodynamic state");
 777         setname[s] = strdup(buf);
 778         s+=1;
 779     }
 780
 781     if (fep->bPrintEnergy)
 782     {
 783         sprintf(buf,"%s (%s)","Energy",unit_energy);
 784         setname[s] = strdup(buf);
 785         s+=1;
 786     }
 787
 788     if (fep->dhdl_derivatives == edhdlderivativesYES)
 789     {
 790         for (i=0;i<efptNR;i++)
 791         {
 792             if (fep->separate_dvdl[i]) {
 793
 794                 if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))
 795                 {
 796                     /* compatibility output */
 797                     sprintf(buf,"%s %s %.4f",dhdl,lambda, fep->init_lambda);
 798                 }
 799                 else
 800                 {
 801                     double lam=fep->init_lambda;
 802                     if (fep->init_lambda < 0)
 803                     {
 804                         lam=fep->all_lambda[i][fep->init_fep_state];
 805                     }
 806                     sprintf(buf,"%s %s = %.4f",dhdl, efpt_singular_names[i],
 807                             lam);
 808                 }
 809                 setname[s] = strdup(buf);
 810                 s+=1;
 811             }
 812         }
 813     }
 814
 815     if (fep->n_lambda > 0)
 816     {
 817         /* g_bar has to determine the lambda values used in this simulation
 818          * from this xvg legend.
 819          */
 820
 821         if (expand->elmcmove > elmcmoveNO) {
 822             nsetsbegin = 1;  /* for including the expanded ensemble */
 823         } else {
 824             nsetsbegin = 0;
 825         }
 826
 827         if (fep->bPrintEnergy)
 828         {
 829             nsetsbegin += 1;
 830         }
 831         nsetsbegin += nsets_dhdl;
 832
 833         for(i=fep->lambda_start_n; i<fep->lambda_stop_n; i++)
 834         {
 835             print_lambda_vector(fep, i, FALSE, FALSE, lambda_vec_str);
 836             if ( (fep->init_lambda >= 0)  && (n_lambda_terms == 1 ))
 837             {
 838                 /* for compatible dhdl.xvg files */
 839                 nps = sprintf(buf,"%s %s %s",deltag,lambda, lambda_vec_str);
 840             }
 841             else
 842             {
 843                 nps = sprintf(buf,"%s %s to %s",deltag,lambda, lambda_vec_str);
 844             }
 845
 846             if (ir->bSimTemp)
 847             {
 848                 /* print the temperature for this state if doing simulated annealing */
 849                 sprintf(&buf[nps],"T = %g (%s)",
 850                         ir->simtempvals->temperatures[s-(nsetsbegin)],
 851                         unit_temp_K);
 852             }
 853             setname[s] = strdup(buf);
 854             s++;
 855         }
 856         if (write_pV) {
 857             np = sprintf(buf,"pV (%s)",unit_energy);
 858             setname[nsetsextend-1] = strdup(buf);  /* the first entry after
 859                                                       nsets */
 860         }
 861
 862         xvgr_legend(fp,nsetsextend,(const char **)setname,oenv);
 863
 864         for(s=0; s<nsetsextend; s++)
 865         {
 866             sfree(setname[s]);
 867         }
 868         sfree(setname);
 869     }
 870
 871     return fp;
 872 }
 873
 874 static void copy_energy(t_mdebin *md, real e[],real ecpy[])
 875 {
 876     int i,j;
 877
 878     for(i=j=0; (i<F_NRE); i++)
 879         if (md->bEner[i])
 880             ecpy[j++] = e[i];
 881     if (j != md->f_nre)
 882         gmx_incons("Number of energy terms wrong");
 883 }
 884
 885 void upd_mdebin(t_mdebin *md,
 886                 gmx_bool bDoDHDL,
 887                 gmx_bool bSum,
 888                 double time,
 889                 real tmass,
 890                 gmx_enerdata_t *enerd,
 891                 t_state *state,
 892                 t_lambda *fep,
 893                 t_expanded *expand,
 894                 matrix  box,
 895                 tensor svir,
 896                 tensor fvir,
 897                 tensor vir,
 898                 tensor pres,
 899                 gmx_ekindata_t *ekind,
 900                 rvec mu_tot,
 901                 gmx_constr_t constr)
 902 {
 903     int    i,j,k,kk,m,n,gid;
 904     real   crmsd[2],tmp6[6];
 905     real   bs[NTRICLBOXS],vol,dens,pv,enthalpy;
 906     real   eee[egNR];
 907     real   ecopy[F_NRE];
 908     double store_dhdl[efptNR];
 909     real   store_energy=0;
 910     real   tmp;
 911
 912     /* Do NOT use the box in the state variable, but the separate box provided
 913      * as an argument. This is because we sometimes need to write the box from
 914      * the last timestep to match the trajectory frames.
 915      */
 916     copy_energy(md, enerd->term,ecopy);
 917     add_ebin(md->ebin,md->ie,md->f_nre,ecopy,bSum);
 918     if (md->nCrmsd)
 919     {
 920         crmsd[0] = constr_rmsd(constr,FALSE);
 921         if (md->nCrmsd > 1)
 922         {
 923             crmsd[1] = constr_rmsd(constr,TRUE);
 924         }
 925         add_ebin(md->ebin,md->iconrmsd,md->nCrmsd,crmsd,FALSE);
 926     }
 927     if (md->bDynBox)
 928     {
 929         int nboxs;
 930         if(md->bTricl)
 931         {
 932             bs[0] = box[XX][XX];
 933             bs[1] = box[YY][YY];
 934             bs[2] = box[ZZ][ZZ];
 935             bs[3] = box[YY][XX];
 936             bs[4] = box[ZZ][XX];
 937             bs[5] = box[ZZ][YY];
 938             nboxs=NTRICLBOXS;
 939         }
 940         else
 941         {
 942             bs[0] = box[XX][XX];
 943             bs[1] = box[YY][YY];
 944             bs[2] = box[ZZ][ZZ];
 945             nboxs=NBOXS;
 946         }
 947         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 948         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 949         add_ebin(md->ebin,md->ib   ,nboxs,bs   ,bSum);
 950         add_ebin(md->ebin,md->ivol ,1    ,&vol ,bSum);
 951         add_ebin(md->ebin,md->idens,1    ,&dens,bSum);
 952
 953         if (md->bDiagPres)
 954         {
 955             /* This is pV (in kJ/mol).  The pressure is the reference pressure,
 956                not the instantaneous pressure */
 957             pv = vol*md->ref_p/PRESFAC;
 958
 959             add_ebin(md->ebin,md->ipv  ,1    ,&pv  ,bSum);
 960             enthalpy = pv + enerd->term[F_ETOT];
 961             add_ebin(md->ebin,md->ienthalpy  ,1    ,&enthalpy  ,bSum);
 962         }
 963     }
 964     if (md->bConstrVir)
 965     {
 966         add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
 967         add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
 968     }
 969     if (md->bVir)
 970         add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
 971     if (md->bPress)
 972         add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
 973     if (md->bSurft){
 974         tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
 975         add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
 976     }
 977     if (md->epc == epcPARRINELLORAHMAN || md->epc == epcMTTK)
 978     {
 979         tmp6[0] = state->boxv[XX][XX];
 980         tmp6[1] = state->boxv[YY][YY];
 981         tmp6[2] = state->boxv[ZZ][ZZ];
 982         tmp6[3] = state->boxv[YY][XX];
 983         tmp6[4] = state->boxv[ZZ][XX];
 984         tmp6[5] = state->boxv[ZZ][YY];
 985         add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
 986     }
 987     if (md->bMu)
 988     {
 989         add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
 990     }
 991     if (ekind && ekind->cosacc.cos_accel != 0)
 992     {
 993         vol  = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
 994         dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
 995         add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
 996         /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
 997         tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
 998                  *dens*vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
 999         add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
1000     }
1001     if (md->nE > 1)
1002     {
1003         n=0;
1004         for(i=0; (i<md->nEg); i++)
1005         {
1006             for(j=i; (j<md->nEg); j++)
1007             {
1008                 gid=GID(i,j,md->nEg);
1009                 for(k=kk=0; (k<egNR); k++)
1010                 {
1011                     if (md->bEInd[k])
1012                     {
1013                         eee[kk++] = enerd->grpp.ener[k][gid];
1014                     }
1015                 }
1016                 add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
1017                 n++;
1018             }
1019         }
1020     }
1021
1022     if (ekind)
1023     {
1024         for(i=0; (i<md->nTC); i++)
1025         {
1026             md->tmp_r[i] = ekind->tcstat[i].T;
1027         }
1028         add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
1029
1030         if (md->etc == etcNOSEHOOVER)
1031         {
1032             /* whether to print Nose-Hoover chains: */
1033             if (md->bPrintNHChains)
1034             {
1035                 if (md->bNHC_trotter)
1036                 {
1037                     for(i=0; (i<md->nTC); i++)
1038                     {
1039                         for (j=0;j<md->nNHC;j++)
1040                         {
1041                             k = i*md->nNHC+j;
1042                             md->tmp_r[2*k] = state->nosehoover_xi[k];
1043                             md->tmp_r[2*k+1] = state->nosehoover_vxi[k];
1044                         }
1045                     }
1046                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
1047
1048                     if (md->bMTTK) {
1049                         for(i=0; (i<md->nTCP); i++)
1050                         {
1051                             for (j=0;j<md->nNHC;j++)
1052                             {
1053                                 k = i*md->nNHC+j;
1054                                 md->tmp_r[2*k] = state->nhpres_xi[k];
1055                                 md->tmp_r[2*k+1] = state->nhpres_vxi[k];
1056                             }
1057                         }
1058                         add_ebin(md->ebin,md->itcb,md->mdeb_n,md->tmp_r,bSum);
1059                     }
1060                 }
1061                 else
1062                 {
1063                     for(i=0; (i<md->nTC); i++)
1064                     {
1065                         md->tmp_r[2*i] = state->nosehoover_xi[i];
1066                         md->tmp_r[2*i+1] = state->nosehoover_vxi[i];
1067                     }
1068                     add_ebin(md->ebin,md->itc,md->mde_n,md->tmp_r,bSum);
1069                 }
1070             }
1071         }
1072         else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
1073                  md->etc == etcVRESCALE)
1074         {
1075             for(i=0; (i<md->nTC); i++)
1076             {
1077                 md->tmp_r[i] = ekind->tcstat[i].lambda;
1078             }
1079             add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
1080         }
1081     }
1082
1083     if (ekind && md->nU > 1)
1084     {
1085         for(i=0; (i<md->nU); i++)
1086         {
1087             copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
1088         }
1089         add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
1090     }
1091
1092     ebin_increase_count(md->ebin,bSum);
1093
1094     /* BAR + thermodynamic integration values */
1095     if ((md->fp_dhdl || md->dhc) && bDoDHDL)
1096     {
1097         for(i=0; i<enerd->n_lambda-1; i++) {
1098             /* zero for simulated tempering */
1099             md->dE[i] = enerd->enerpart_lambda[i+1]-enerd->enerpart_lambda[0];
1100             if (md->temperatures!=NULL)
1101             {
1102                 /* MRS: is this right, given the way we have defined the exchange probabilities? */
1103                 /* is this even useful to have at all? */
1104                 md->dE[i] += (md->temperatures[i]/
1105                           md->temperatures[state->fep_state]-1.0)*
1106                             enerd->term[F_EKIN];
1107             }
1108         }
1109
1110         if (md->fp_dhdl)
1111         {
1112             fprintf(md->fp_dhdl,"%.4f",time);
1113             /* the current free energy state */
1114
1115             /* print the current state if we are doing expanded ensemble */
1116             if (expand->elmcmove > elmcmoveNO) {
1117                 fprintf(md->fp_dhdl," %4d",state->fep_state);
1118             }
1119             /* total energy (for if the temperature changes */
1120             if (fep->bPrintEnergy)
1121             {
1122                 store_energy = enerd->term[F_ETOT];
1123                 fprintf(md->fp_dhdl," %#.8g",store_energy);
1124             }
1125
1126             if (fep->dhdl_derivatives == edhdlderivativesYES)
1127             {
1128                 for (i=0;i<efptNR;i++)
1129                 {
1130                     if (fep->separate_dvdl[i])
1131                     {
1132                         /* assumes F_DVDL is first */
1133                         fprintf(md->fp_dhdl," %#.8g",enerd->term[F_DVDL+i]);
1134                     }
1135                 }
1136             }
1137             for(i=fep->lambda_start_n;i<fep->lambda_stop_n;i++)
1138             {
1139                 fprintf(md->fp_dhdl," %#.8g",md->dE[i]);
1140             }
1141             if ((md->epc!=epcNO)  &&
1142                 (enerd->n_lambda > 0) &&
1143                 (fep->init_lambda<0))
1144             {
1145                 fprintf(md->fp_dhdl," %#.8g",pv);  /* PV term only needed when
1146                                                       there are alternate state
1147                                                       lambda and we're not in
1148                                                       compatibility mode */
1149             }
1150             fprintf(md->fp_dhdl,"\n");
1151             /* and the binary free energy output */
1152         }
1153         if (md->dhc && bDoDHDL)
1154         {
1155             int idhdl = 0;
1156             for (i=0;i<efptNR;i++)
1157             {
1158                 if (fep->separate_dvdl[i])
1159                 {
1160                     /* assumes F_DVDL is first */
1161                     store_dhdl[idhdl] = enerd->term[F_DVDL+i];
1162                     idhdl+=1;
1163                 }
1164             }
1165             store_energy = enerd->term[F_ETOT];
1166             /* store_dh is dE */
1167             mde_delta_h_coll_add_dh(md->dhc,
1168                                     (double)state->fep_state,
1169                                     store_energy,
1170                                     pv,
1171                                     store_dhdl,
1172                                     md->dE + fep->lambda_start_n,
1173                                     time);
1174         }
1175     }
1176 }
1177
1178
1179 void upd_mdebin_step(t_mdebin *md)
1180 {
1181     ebin_increase_count(md->ebin,FALSE);
1182 }
1183
1184 static void npr(FILE *log,int n,char c)
1185 {
1186     for(; (n>0); n--) fprintf(log,"%c",c);
1187 }
1188
1189 static void pprint(FILE *log,const char *s,t_mdebin *md)
1190 {
1191     char CHAR='#';
1192     int  slen;
1193     char buf1[22],buf2[22];
1194
1195     slen = strlen(s);
1196     fprintf(log,"\t<======  ");
1197     npr(log,slen,CHAR);
1198     fprintf(log,"  ==>\n");
1199     fprintf(log,"\t<====  %s  ====>\n",s);
1200     fprintf(log,"\t<==  ");
1201     npr(log,slen,CHAR);
1202     fprintf(log,"  ======>\n\n");
1203
1204     fprintf(log,"\tStatistics over %s steps using %s frames\n",
1205             gmx_step_str(md->ebin->nsteps_sim,buf1),
1206             gmx_step_str(md->ebin->nsum_sim,buf2));
1207     fprintf(log,"\n");
1208 }
1209
1210 void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lambda)
1211 {
1212     char buf[22];
1213
1214     fprintf(log,"   %12s   %12s   %12s\n"
1215             "   %12s   %12.5f   %12.5f\n\n",
1216             "Step","Time","Lambda",gmx_step_str(steps,buf),time,lambda);
1217 }
1218
1219 void print_ebin(ener_file_t fp_ene,gmx_bool bEne,gmx_bool bDR,gmx_bool bOR,
1220                 FILE *log,
1221                 gmx_large_int_t step,double time,
1222                 int mode,gmx_bool bCompact,
1223                 t_mdebin *md,t_fcdata *fcd,
1224                 gmx_groups_t *groups,t_grpopts *opts)
1225 {
1226     /*static char **grpnms=NULL;*/
1227     char        buf[246];
1228     int         i,j,n,ni,nj,ndr,nor,b;
1229     int         ndisre=0;
1230     real        *disre_rm3tav, *disre_rt;
1231
1232     /* these are for the old-style blocks (1 subblock, only reals), because
1233        there can be only one per ID for these */
1234     int         nr[enxNR];
1235     int         id[enxNR];
1236     real        *block[enxNR];
1237
1238     /* temporary arrays for the lambda values to write out */
1239     double      enxlambda_data[2];
1240
1241     t_enxframe  fr;
1242
1243     switch (mode)
1244     {
1245         case eprNORMAL:
1246             init_enxframe(&fr);
1247             fr.t            = time;
1248             fr.step         = step;
1249             fr.nsteps       = md->ebin->nsteps;
1250             fr.dt           = md->delta_t;
1251             fr.nsum         = md->ebin->nsum;
1252             fr.nre          = (bEne) ? md->ebin->nener : 0;
1253             fr.ener         = md->ebin->e;
1254             ndisre          = bDR ? fcd->disres.npair : 0;
1255             disre_rm3tav    = fcd->disres.rm3tav;
1256             disre_rt        = fcd->disres.rt;
1257             /* Optional additional old-style (real-only) blocks. */
1258             for(i=0; i<enxNR; i++)
1259             {
1260                 nr[i] = 0;
1261             }
1262             if (fcd->orires.nr > 0 && bOR)
1263             {
1264                 diagonalize_orires_tensors(&(fcd->orires));
1265                 nr[enxOR]     = fcd->orires.nr;
1266                 block[enxOR]  = fcd->orires.otav;
1267                 id[enxOR]     = enxOR;
1268                 nr[enxORI]    = (fcd->orires.oinsl != fcd->orires.otav) ?
1269                           fcd->orires.nr : 0;
1270                 block[enxORI] = fcd->orires.oinsl;
1271                 id[enxORI]    = enxORI;
1272                 nr[enxORT]    = fcd->orires.nex*12;
1273                 block[enxORT] = fcd->orires.eig;
1274                 id[enxORT]    = enxORT;
1275             }
1276
1277             /* whether we are going to wrte anything out: */
1278             if (fr.nre || ndisre || nr[enxOR] || nr[enxORI])
1279             {
1280
1281                 /* the old-style blocks go first */
1282                 fr.nblock = 0;
1283                 for(i=0; i<enxNR; i++)
1284                 {
1285                     if (nr[i] > 0)
1286                     {
1287                         fr.nblock = i + 1;
1288                     }
1289                 }
1290                 add_blocks_enxframe(&fr, fr.nblock);
1291                 for(b=0;b<fr.nblock;b++)
1292                 {
1293                     add_subblocks_enxblock(&(fr.block[b]), 1);
1294                     fr.block[b].id=id[b];
1295                     fr.block[b].sub[0].nr = nr[b];
1296 #ifndef GMX_DOUBLE
1297                     fr.block[b].sub[0].type = xdr_datatype_float;
1298                     fr.block[b].sub[0].fval = block[b];
1299 #else
1300                     fr.block[b].sub[0].type = xdr_datatype_double;
1301                     fr.block[b].sub[0].dval = block[b];
1302 #endif
1303                 }
1304
1305                 /* check for disre block & fill it. */
1306                 if (ndisre>0)
1307                 {
1308                     int db = fr.nblock;
1309                     fr.nblock+=1;
1310                     add_blocks_enxframe(&fr, fr.nblock);
1311
1312                     add_subblocks_enxblock(&(fr.block[db]), 2);
1313                     fr.block[db].id=enxDISRE;
1314                     fr.block[db].sub[0].nr=ndisre;
1315                     fr.block[db].sub[1].nr=ndisre;
1316 #ifndef GMX_DOUBLE
1317                     fr.block[db].sub[0].type=xdr_datatype_float;
1318                     fr.block[db].sub[1].type=xdr_datatype_float;
1319                     fr.block[db].sub[0].fval=disre_rt;
1320                     fr.block[db].sub[1].fval=disre_rm3tav;
1321 #else
1322                     fr.block[db].sub[0].type=xdr_datatype_double;
1323                     fr.block[db].sub[1].type=xdr_datatype_double;
1324                     fr.block[db].sub[0].dval=disre_rt;
1325                     fr.block[db].sub[1].dval=disre_rm3tav;
1326 #endif
1327                 }
1328                 /* here we can put new-style blocks */
1329
1330                 /* Free energy perturbation blocks */
1331                 if (md->dhc)
1332                 {
1333                     mde_delta_h_coll_handle_block(md->dhc, &fr, fr.nblock);
1334                 }
1335
1336                 /* we can now free & reset the data in the blocks */
1337                 if (md->dhc)
1338                 {
1339                     mde_delta_h_coll_reset(md->dhc);
1340                 }
1341
1342                 /* do the actual I/O */
1343                 do_enx(fp_ene,&fr);
1344                 gmx_fio_check_file_position(enx_file_pointer(fp_ene));
1345                 if (fr.nre)
1346                 {
1347                     /* We have stored the sums, so reset the sum history */
1348                     reset_ebin_sums(md->ebin);
1349                 }
1350             }
1351             free_enxframe(&fr);
1352             break;
1353         case eprAVER:
1354             if (log)
1355             {
1356                 pprint(log,"A V E R A G E S",md);
1357             }
1358             break;
1359         case eprRMS:
1360             if (log)
1361             {
1362                 pprint(log,"R M S - F L U C T U A T I O N S",md);
1363             }
1364             break;
1365         default:
1366             gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
1367     }
1368
1369     if (log)
1370     {
1371         for(i=0;i<opts->ngtc;i++)
1372         {
1373             if(opts->annealing[i]!=eannNO)
1374             {
1375                 fprintf(log,"Current ref_t for group %s: %8.1f\n",
1376                         *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
1377                         opts->ref_t[i]);
1378             }
1379         }
1380         if (mode==eprNORMAL && fcd->orires.nr>0)
1381         {
1382             print_orires_log(log,&(fcd->orires));
1383         }
1384         fprintf(log,"   Energies (%s)\n",unit_energy);
1385         pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
1386         fprintf(log,"\n");
1387
1388         if (!bCompact)
1389         {
1390             if (md->bDynBox)
1391             {
1392                 pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
1393                         mode,TRUE);
1394                 fprintf(log,"\n");
1395             }
1396             if (md->bConstrVir)
1397             {
1398                 fprintf(log,"   Constraint Virial (%s)\n",unit_energy);
1399                 pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
1400                 fprintf(log,"\n");
1401                 fprintf(log,"   Force Virial (%s)\n",unit_energy);
1402                 pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
1403                 fprintf(log,"\n");
1404             }
1405             if (md->bVir)
1406             {
1407                 fprintf(log,"   Total Virial (%s)\n",unit_energy);
1408                 pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
1409                 fprintf(log,"\n");
1410             }
1411             if (md->bPress)
1412             {
1413                 fprintf(log,"   Pressure (%s)\n",unit_pres_bar);
1414                 pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
1415                 fprintf(log,"\n");
1416             }
1417             if (md->bMu)
1418             {
1419                 fprintf(log,"   Total Dipole (%s)\n",unit_dipole_D);
1420                 pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
1421                 fprintf(log,"\n");
1422             }
1423
1424             if (md->nE > 1)
1425             {
1426                 if (md->print_grpnms==NULL)
1427                 {
1428                     snew(md->print_grpnms,md->nE);
1429                     n=0;
1430                     for(i=0; (i<md->nEg); i++)
1431                     {
1432                         ni=groups->grps[egcENER].nm_ind[i];
1433                         for(j=i; (j<md->nEg); j++)
1434                         {
1435                             nj=groups->grps[egcENER].nm_ind[j];
1436                             sprintf(buf,"%s-%s",*(groups->grpname[ni]),
1437                                     *(groups->grpname[nj]));
1438                             md->print_grpnms[n++]=strdup(buf);
1439                         }
1440                     }
1441                 }
1442                 sprintf(buf,"Epot (%s)",unit_energy);
1443                 fprintf(log,"%15s   ",buf);
1444                 for(i=0; (i<egNR); i++)
1445                 {
1446                     if (md->bEInd[i])
1447                     {
1448                         fprintf(log,"%12s   ",egrp_nm[i]);
1449                     }
1450                 }
1451                 fprintf(log,"\n");
1452                 for(i=0; (i<md->nE); i++)
1453                 {
1454                     fprintf(log,"%15s",md->print_grpnms[i]);
1455                     pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
1456                             FALSE);
1457                 }
1458                 fprintf(log,"\n");
1459             }
1460             if (md->nTC > 1)
1461             {
1462                 pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
1463                 fprintf(log,"\n");
1464             }
1465             if (md->nU > 1)
1466             {
1467                 fprintf(log,"%15s   %12s   %12s   %12s\n",
1468                         "Group","Ux","Uy","Uz");
1469                 for(i=0; (i<md->nU); i++)
1470                 {
1471                     ni=groups->grps[egcACC].nm_ind[i];
1472                     fprintf(log,"%15s",*groups->grpname[ni]);
1473                     pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
1474                 }
1475                 fprintf(log,"\n");
1476             }
1477         }
1478     }
1479
1480 }
1481
1482 void update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
1483 {
1484     int i;
1485
1486     enerhist->nsteps     = mdebin->ebin->nsteps;
1487     enerhist->nsum       = mdebin->ebin->nsum;
1488     enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
1489     enerhist->nsum_sim   = mdebin->ebin->nsum_sim;
1490     enerhist->nener      = mdebin->ebin->nener;
1491
1492     if (mdebin->ebin->nsum > 0)
1493     {
1494         /* Check if we need to allocate first */
1495         if(enerhist->ener_ave == NULL)
1496         {
1497             snew(enerhist->ener_ave,enerhist->nener);
1498             snew(enerhist->ener_sum,enerhist->nener);
1499         }
1500
1501         for(i=0;i<enerhist->nener;i++)
1502         {
1503             enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
1504             enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
1505         }
1506     }
1507
1508     if (mdebin->ebin->nsum_sim > 0)
1509     {
1510         /* Check if we need to allocate first */
1511         if(enerhist->ener_sum_sim == NULL)
1512         {
1513             snew(enerhist->ener_sum_sim,enerhist->nener);
1514         }
1515
1516         for(i=0;i<enerhist->nener;i++)
1517         {
1518             enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
1519         }
1520     }
1521     if (mdebin->dhc)
1522     {
1523         mde_delta_h_coll_update_energyhistory(mdebin->dhc, enerhist);
1524     }
1525 }
1526
1527 void restore_energyhistory_from_state(t_mdebin * mdebin,
1528                                       energyhistory_t * enerhist)
1529 {
1530     int i;
1531
1532     if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
1533         mdebin->ebin->nener != enerhist->nener)
1534     {
1535         gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
1536                   mdebin->ebin->nener,enerhist->nener);
1537     }
1538
1539     mdebin->ebin->nsteps     = enerhist->nsteps;
1540     mdebin->ebin->nsum       = enerhist->nsum;
1541     mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
1542     mdebin->ebin->nsum_sim   = enerhist->nsum_sim;
1543
1544     for(i=0; i<mdebin->ebin->nener; i++)
1545     {
1546         mdebin->ebin->e[i].eav  =
1547                   (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
1548         mdebin->ebin->e[i].esum =
1549                   (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
1550         mdebin->ebin->e_sim[i].esum =
1551                   (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
1552     }
1553     if (mdebin->dhc)
1554     {
1555         mde_delta_h_coll_restore_energyhistory(mdebin->dhc, enerhist);
1556     }
1557 }