src/mdlib/mdebin.c

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