src/tools/gmx_energy.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  * Green Red Orange Magenta Azure Cyan Skyblue
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <string.h>
  41 #include <math.h>
  42 #include <ctype.h>
  43
  44 #include "typedefs.h"
  45 #include "gmx_fatal.h"
  46 #include "vec.h"
  47 #include "string2.h"
  48 #include "smalloc.h"
  49 #include "enxio.h"
  50 #include "statutil.h"
  51 #include "names.h"
  52 #include "copyrite.h"
  53 #include "macros.h"
  54 #include "xvgr.h"
  55 #include "gstat.h"
  56 #include "physics.h"
  57 #include "tpxio.h"
  58 #include "viewit.h"
  59 #include "mtop_util.h"
  60 #include "gmx_ana.h"
  61 #include "mdebin.h"
  62
  63 static real       minthird=-1.0/3.0,minsixth=-1.0/6.0;
  64
  65 typedef struct {
  66     real sum;
  67     real sum2;
  68 } exactsum_t;
  69
  70 typedef struct {
  71     real   *ener;
  72     exactsum_t *es;
  73     gmx_bool   bExactStat;
  74     double av;
  75     double rmsd;
  76     double ee;
  77     double slope;
  78 } enerdat_t;
  79
  80 typedef struct {
  81     gmx_large_int_t nsteps;
  82     gmx_large_int_t npoints;
  83     int             nframes;
  84     int             *step;
  85     int             *steps;
  86     int             *points;
  87     enerdat_t       *s;
  88 } enerdata_t;
  89
  90 static double mypow(double x,double y)
  91 {
  92   if (x > 0)
  93     return pow(x,y);
  94   else
  95     return 0.0;
  96 }
  97
  98 static int *select_it(int nre,char *nm[],int *nset)
  99 {
 100   gmx_bool *bE;
 101   int  n,k,j,i;
 102   int  *set;
 103   gmx_bool bVerbose = TRUE;
 104
 105   if ((getenv("VERBOSE")) != NULL)
 106     bVerbose = FALSE;
 107
 108   fprintf(stderr,"Select the terms you want from the following list\n");
 109   fprintf(stderr,"End your selection with 0\n");
 110
 111   if ( bVerbose ) {
 112     for(k=0; (k<nre); ) {
 113       for(j=0; (j<4) && (k<nre); j++,k++)
 114         fprintf(stderr," %3d=%14s",k+1,nm[k]);
 115       fprintf(stderr,"\n");
 116     }
 117   }
 118
 119   snew(bE,nre);
 120   do {
 121     if(1 != scanf("%d",&n))
 122     {
 123       gmx_fatal(FARGS,"Error reading user input");
 124     }
 125     if ((n>0) && (n<=nre))
 126       bE[n-1]=TRUE;
 127   } while (n != 0);
 128
 129   snew(set,nre);
 130   for(i=(*nset)=0; (i<nre); i++)
 131     if (bE[i])
 132       set[(*nset)++]=i;
 133
 134   sfree(bE);
 135
 136   return set;
 137 }
 138
 139 static void chomp(char *buf)
 140 {
 141   int len = strlen(buf);
 142
 143   while ((len > 0) && (buf[len-1] == '\n')) {
 144     buf[len-1] = '\0';
 145     len--;
 146   }
 147 }
 148
 149 static int *select_by_name(int nre,gmx_enxnm_t *nm,int *nset)
 150 {
 151   gmx_bool *bE;
 152   int  n,k,kk,j,i,nmatch,nind,nss;
 153   int  *set;
 154   gmx_bool bEOF,bVerbose = TRUE,bLong=FALSE;
 155   char *ptr,buf[STRLEN];
 156   const char *fm4="%3d  %-14s";
 157   const char *fm2="%3d  %-34s";
 158   char **newnm=NULL;
 159
 160   if ((getenv("VERBOSE")) != NULL)
 161     bVerbose = FALSE;
 162
 163   fprintf(stderr,"\n");
 164   fprintf(stderr,"Select the terms you want from the following list by\n");
 165   fprintf(stderr,"selecting either (part of) the name or the number or a combination.\n");
 166   fprintf(stderr,"End your selection with an empty line or a zero.\n");
 167   fprintf(stderr,"-------------------------------------------------------------------\n");
 168
 169   snew(newnm,nre);
 170   j = 0;
 171   for(k=0; k<nre; k++) {
 172     newnm[k] = strdup(nm[k].name);
 173     /* Insert dashes in all the names */
 174     while ((ptr = strchr(newnm[k],' ')) != NULL) {
 175       *ptr = '-';
 176     }
 177     if ( bVerbose ) {
 178       if (j == 0) {
 179         if (k > 0) {
 180           fprintf(stderr,"\n");
 181         }
 182         bLong = FALSE;
 183         for(kk=k; kk<k+4; kk++) {
 184           if (kk < nre && strlen(nm[kk].name) > 14) {
 185             bLong = TRUE;
 186           }
 187         }
 188       } else {
 189         fprintf(stderr," ");
 190       }
 191       if (!bLong) {
 192         fprintf(stderr,fm4,k+1,newnm[k]);
 193         j++;
 194         if (j == 4) {
 195           j = 0;
 196         }
 197       } else {
 198         fprintf(stderr,fm2,k+1,newnm[k]);
 199         j++;
 200         if (j == 2) {
 201           j = 0;
 202         }
 203       }
 204     }
 205   }
 206   if ( bVerbose ) {
 207     fprintf(stderr,"\n\n");
 208   }
 209
 210   snew(bE,nre);
 211
 212   bEOF = FALSE;
 213   while (!bEOF && (fgets2(buf,STRLEN-1,stdin))) {
 214     /* Remove newlines */
 215     chomp(buf);
 216
 217     /* Remove spaces */
 218     trim(buf);
 219
 220     /* Empty line means end of input */
 221     bEOF = (strlen(buf) == 0);
 222     if (!bEOF) {
 223       ptr = buf;
 224       do {
 225         if (!bEOF) {
 226           /* First try to read an integer */
 227           nss   = sscanf(ptr,"%d",&nind);
 228           if (nss == 1) {
 229             /* Zero means end of input */
 230             if (nind == 0) {
 231               bEOF = TRUE;
 232             } else if ((1<=nind) && (nind<=nre)) {
 233               bE[nind-1] = TRUE;
 234             } else {
 235               fprintf(stderr,"number %d is out of range\n",nind);
 236             }
 237           }
 238           else {
 239             /* Now try to read a string */
 240             i = strlen(ptr);
 241             nmatch = 0;
 242             for(nind=0; nind<nre; nind++) {
 243               if (gmx_strcasecmp(newnm[nind],ptr) == 0) {
 244                 bE[nind] = TRUE;
 245                 nmatch++;
 246               }
 247             }
 248             if (nmatch == 0) {
 249               i = strlen(ptr);
 250               nmatch = 0;
 251               for(nind=0; nind<nre; nind++) {
 252                 if (gmx_strncasecmp(newnm[nind],ptr,i) == 0) {
 253                   bE[nind] = TRUE;
 254                   nmatch++;
 255                 }
 256               }
 257               if (nmatch == 0) {
 258                 fprintf(stderr,"String '%s' does not match anything\n",ptr);
 259               }
 260             }
 261           }
 262         }
 263         /* Look for the first space, and remove spaces from there */
 264         if ((ptr = strchr(ptr,' ')) != NULL) {
 265           trim(ptr);
 266         }
 267       } while (!bEOF && (ptr && (strlen(ptr) > 0)));
 268     }
 269   }
 270
 271   snew(set,nre);
 272   for(i=(*nset)=0; (i<nre); i++)
 273     if (bE[i])
 274       set[(*nset)++]=i;
 275
 276   sfree(bE);
 277
 278   if (*nset == 0)
 279     gmx_fatal(FARGS,"No energy terms selected");
 280
 281   for(i=0; (i<nre); i++)
 282     sfree(newnm[i]);
 283   sfree(newnm);
 284
 285   return set;
 286 }
 287
 288 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
 289 {
 290     gmx_mtop_t mtop;
 291     int        natoms;
 292     t_iatom    *iatom;
 293     matrix     box;
 294
 295     /* all we need is the ir to be able to write the label */
 296     read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,&mtop);
 297 }
 298
 299 static void get_orires_parms(const char *topnm,
 300                                int *nor,int *nex,int **label,real **obs)
 301 {
 302   gmx_mtop_t mtop;
 303   gmx_localtop_t *top;
 304   t_inputrec ir;
 305   t_iparams  *ip;
 306   int        natoms,i;
 307   t_iatom    *iatom;
 308   int        nb;
 309   matrix     box;
 310
 311   read_tpx(topnm,&ir,box,&natoms,NULL,NULL,NULL,&mtop);
 312   top = gmx_mtop_generate_local_top(&mtop,&ir);
 313
 314   ip       = top->idef.iparams;
 315   iatom    = top->idef.il[F_ORIRES].iatoms;
 316
 317   /* Count how many distance restraint there are... */
 318   nb = top->idef.il[F_ORIRES].nr;
 319   if (nb == 0)
 320     gmx_fatal(FARGS,"No orientation restraints in topology!\n");
 321
 322   *nor = nb/3;
 323   *nex = 0;
 324   snew(*label,*nor);
 325   snew(*obs,*nor);
 326   for(i=0; i<nb; i+=3) {
 327     (*label)[i/3] = ip[iatom[i]].orires.label;
 328     (*obs)[i/3]   = ip[iatom[i]].orires.obs;
 329     if (ip[iatom[i]].orires.ex >= *nex)
 330       *nex = ip[iatom[i]].orires.ex+1;
 331   }
 332   fprintf(stderr,"Found %d orientation restraints with %d experiments",
 333           *nor,*nex);
 334 }
 335
 336 static int get_bounds(const char *topnm,
 337                       real **bounds,int **index,int **dr_pair,int *npairs,
 338                       gmx_mtop_t *mtop,gmx_localtop_t **ltop,t_inputrec *ir)
 339 {
 340   gmx_localtop_t *top;
 341   t_functype *functype;
 342   t_iparams  *ip;
 343   int        natoms,i,j,k,type,ftype,natom;
 344   t_ilist    *disres;
 345   t_iatom    *iatom;
 346   real       *b;
 347   int        *ind,*pair;
 348   int        nb,label1;
 349   matrix     box;
 350
 351   read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,mtop);
 352   snew(*ltop,1);
 353   top = gmx_mtop_generate_local_top(mtop,ir);
 354   *ltop = top;
 355
 356   functype = top->idef.functype;
 357   ip       = top->idef.iparams;
 358
 359   /* Count how many distance restraint there are... */
 360   nb=top->idef.il[F_DISRES].nr;
 361   if (nb == 0)
 362     gmx_fatal(FARGS,"No distance restraints in topology!\n");
 363
 364   /* Allocate memory */
 365   snew(b,nb);
 366   snew(ind,nb);
 367   snew(pair,nb+1);
 368
 369   /* Fill the bound array */
 370   nb=0;
 371   for(i=0; (i<top->idef.ntypes); i++) {
 372     ftype = functype[i];
 373     if (ftype == F_DISRES) {
 374
 375       label1 = ip[i].disres.label;
 376       b[nb]   = ip[i].disres.up1;
 377       ind[nb] = label1;
 378       nb++;
 379     }
 380   }
 381   *bounds = b;
 382
 383   /* Fill the index array */
 384   label1  = -1;
 385   disres  = &(top->idef.il[F_DISRES]);
 386   iatom   = disres->iatoms;
 387   for(i=j=k=0; (i<disres->nr); ) {
 388     type  = iatom[i];
 389     ftype = top->idef.functype[type];
 390     natom = interaction_function[ftype].nratoms+1;
 391     if (label1 != top->idef.iparams[type].disres.label) {
 392       pair[j] = k;
 393       label1  = top->idef.iparams[type].disres.label;
 394       j ++;
 395     }
 396     k++;
 397     i += natom;
 398   }
 399   pair[j]  = k;
 400   *npairs = k;
 401   if (j != nb)
 402     gmx_incons("get_bounds for distance restraints");
 403
 404   *index   = ind;
 405   *dr_pair = pair;
 406
 407   return nb;
 408 }
 409
 410 static void calc_violations(real rt[],real rav3[],int nb,int index[],
 411                             real bounds[],real *viol,double *st,double *sa)
 412 {
 413   const   real sixth=1.0/6.0;
 414   int     i,j;
 415   double  rsum,rav,sumaver,sumt;
 416
 417   sumaver = 0;
 418   sumt    = 0;
 419   for(i=0; (i<nb); i++) {
 420     rsum = 0.0;
 421     rav  = 0.0;
 422     for(j=index[i]; (j<index[i+1]); j++) {
 423       if (viol)
 424         viol[j] += mypow(rt[j],-3.0);
 425       rav     += sqr(rav3[j]);
 426       rsum    += mypow(rt[j],-6);
 427     }
 428     rsum    = max(0.0,mypow(rsum,-sixth)-bounds[i]);
 429     rav     = max(0.0,mypow(rav, -sixth)-bounds[i]);
 430
 431     sumt    += rsum;
 432     sumaver += rav;
 433   }
 434   *st = sumt;
 435   *sa = sumaver;
 436 }
 437
 438 static void analyse_disre(const char *voutfn,    int nframes,
 439                           real violaver[], real bounds[], int index[],
 440                           int pair[],      int nbounds,
 441                           const output_env_t oenv)
 442 {
 443     FILE   *vout;
 444     double sum,sumt,sumaver;
 445     int    i,j;
 446
 447     /* Subtract bounds from distances, to calculate violations */
 448     calc_violations(violaver,violaver,
 449             nbounds,pair,bounds,NULL,&sumt,&sumaver);
 450
 451 #ifdef DEBUG
 452     fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
 453             sumaver);
 454     fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",
 455             sumt);
 456 #endif
 457     vout=xvgropen(voutfn,"r\\S-3\\N average violations","DR Index","nm",
 458             oenv);
 459     sum  = 0.0;
 460     sumt = 0.0;
 461     for(i=0; (i<nbounds); i++) {
 462         /* Do ensemble averaging */
 463         sumaver = 0;
 464         for(j=pair[i]; (j<pair[i+1]); j++)
 465             sumaver += sqr(violaver[j]/nframes);
 466         sumaver = max(0.0,mypow(sumaver,minsixth)-bounds[i]);
 467
 468         sumt   += sumaver;
 469         sum     = max(sum,sumaver);
 470         fprintf(vout,"%10d  %10.5e\n",index[i],sumaver);
 471     }
 472 #ifdef DEBUG
 473     for(j=0; (j<dr.ndr); j++)
 474         fprintf(vout,"%10d  %10.5e\n",j,mypow(violaver[j]/nframes,minthird));
 475 #endif
 476     ffclose(vout);
 477
 478     fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
 479             sumt);
 480     fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",sum);
 481
 482     do_view(oenv,voutfn,"-graphtype bar");
 483 }
 484
 485 static void einstein_visco(const char *fn,const char *fni,int nsets,
 486                            int nframes,real **sum,
 487                            real V,real T,int nsteps,double time[],
 488                            const output_env_t oenv)
 489 {
 490     FILE *fp0,*fp1;
 491     double av[4],avold[4];
 492     double fac,dt,di;
 493     int  i,j,m,nf4;
 494
 495     if (nframes < 1)
 496         return;
 497
 498     dt  = (time[1]-time[0]);
 499     nf4 = nframes/4+1;
 500
 501     for(i=0; i<=nsets; i++)
 502         avold[i] = 0;
 503     fp0=xvgropen(fni,"Shear viscosity integral",
 504             "Time (ps)","(kg m\\S-1\\N s\\S-1\\N ps)",oenv);
 505     fp1=xvgropen(fn,"Shear viscosity using Einstein relation",
 506             "Time (ps)","(kg m\\S-1\\N s\\S-1\\N)",oenv);
 507     for(i=1; i<nf4; i++) {
 508         fac = dt*nframes/nsteps;
 509         for(m=0; m<=nsets; m++)
 510             av[m] = 0;
 511         for(j=0; j<nframes-i; j++) {
 512             for(m=0; m<nsets; m++) {
 513                 di   = sqr(fac*(sum[m][j+i]-sum[m][j]));
 514
 515                 av[m]     += di;
 516                 av[nsets] += di/nsets;
 517             }
 518         }
 519         /* Convert to SI for the viscosity */
 520         fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
 521         fprintf(fp0,"%10g",time[i]-time[0]);
 522         for(m=0; (m<=nsets); m++) {
 523             av[m] = fac*av[m];
 524             fprintf(fp0,"  %10g",av[m]);
 525         }
 526         fprintf(fp0,"\n");
 527         fprintf(fp1,"%10g",0.5*(time[i]+time[i-1])-time[0]);
 528         for(m=0; (m<=nsets); m++) {
 529             fprintf(fp1,"  %10g",(av[m]-avold[m])/dt);
 530             avold[m] = av[m];
 531         }
 532         fprintf(fp1,"\n");
 533     }
 534     ffclose(fp0);
 535     ffclose(fp1);
 536 }
 537
 538 typedef struct {
 539     gmx_large_int_t np;
 540     double sum;
 541     double sav;
 542     double sav2;
 543 } ee_sum_t;
 544
 545 typedef struct {
 546     int b;
 547     ee_sum_t sum;
 548     gmx_large_int_t nst;
 549     gmx_large_int_t nst_min;
 550 } ener_ee_t;
 551
 552 static void clear_ee_sum(ee_sum_t *ees)
 553 {
 554     ees->sav  = 0;
 555     ees->sav2 = 0;
 556     ees->np   = 0;
 557     ees->sum  = 0;
 558 }
 559
 560 static void add_ee_sum(ee_sum_t *ees,double sum,int np)
 561 {
 562     ees->np  += np;
 563     ees->sum += sum;
 564 }
 565
 566 static void add_ee_av(ee_sum_t *ees)
 567 {
 568     double av;
 569
 570     av = ees->sum/ees->np;
 571     ees->sav  += av;
 572     ees->sav2 += av*av;
 573     ees->np    = 0;
 574     ees->sum   = 0;
 575 }
 576
 577 static double calc_ee2(int nb,ee_sum_t *ees)
 578 {
 579     return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
 580 }
 581
 582 static void set_ee_av(ener_ee_t *eee)
 583 {
 584     if (debug)
 585     {
 586         char buf[STEPSTRSIZE];
 587         fprintf(debug,"Storing average for err.est.: %s steps\n",
 588                 gmx_step_str(eee->nst,buf));
 589     }
 590     add_ee_av(&eee->sum);
 591     eee->b++;
 592     if (eee->b == 1 || eee->nst < eee->nst_min)
 593     {
 594         eee->nst_min = eee->nst;
 595     }
 596     eee->nst = 0;
 597 }
 598
 599 static void calc_averages(int nset,enerdata_t *edat,int nbmin,int nbmax)
 600 {
 601     int  nb,i,f,nee;
 602     double sum,sum2,sump,see2;
 603     gmx_large_int_t steps,np,p,bound_nb;
 604     enerdat_t *ed;
 605     exactsum_t *es;
 606     gmx_bool bAllZero;
 607     double x,sx,sy,sxx,sxy;
 608     ener_ee_t *eee;
 609
 610     /* Check if we have exact statistics over all points */
 611     for(i=0; i<nset; i++)
 612     {
 613         ed = &edat->s[i];
 614         ed->bExactStat = FALSE;
 615         if (edat->npoints > 0)
 616         {
 617             /* All energy file sum entries 0 signals no exact sums.
 618              * But if all energy values are 0, we still have exact sums.
 619              */
 620             bAllZero = TRUE;
 621             for(f=0; f<edat->nframes && !ed->bExactStat; f++)
 622             {
 623                 if (ed->ener[i] != 0)
 624                 {
 625                     bAllZero = FALSE;
 626                 }
 627                 ed->bExactStat = (ed->es[f].sum != 0);
 628             }
 629             if (bAllZero)
 630             {
 631                 ed->bExactStat = TRUE;
 632             }
 633         }
 634     }
 635
 636     snew(eee,nbmax+1);
 637     for(i=0; i<nset; i++)
 638     {
 639         ed = &edat->s[i];
 640
 641         sum  = 0;
 642         sum2 = 0;
 643         np   = 0;
 644         sx   = 0;
 645         sy   = 0;
 646         sxx  = 0;
 647         sxy  = 0;
 648         for(nb=nbmin; nb<=nbmax; nb++)
 649         {
 650             eee[nb].b     = 0;
 651             clear_ee_sum(&eee[nb].sum);
 652             eee[nb].nst   = 0;
 653             eee[nb].nst_min = 0;
 654         }
 655         for(f=0; f<edat->nframes; f++)
 656         {
 657             es = &ed->es[f];
 658
 659             if (ed->bExactStat)
 660             {
 661                 /* Add the sum and the sum of variances to the totals. */
 662                 p     = edat->points[f];
 663                 sump  = es->sum;
 664                 sum2 += es->sum2;
 665                 if (np > 0)
 666                 {
 667                     sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
 668                         *np*(np + p)/p;
 669                 }
 670             }
 671             else
 672             {
 673                 /* Add a single value to the sum and sum of squares. */
 674                 p     = 1;
 675                 sump  = ed->ener[f];
 676                 sum2 += dsqr(sump);
 677             }
 678
 679             /* sum has to be increased after sum2 */
 680             np  += p;
 681             sum += sump;
 682
 683             /* For the linear regression use variance 1/p.
 684              * Note that sump is the sum, not the average, so we don't need p*.
 685              */
 686             x    = edat->step[f] - 0.5*(edat->steps[f] - 1);
 687             sx  += p*x;
 688             sy  += sump;
 689             sxx += p*x*x;
 690             sxy += x*sump;
 691
 692             for(nb=nbmin; nb<=nbmax; nb++)
 693             {
 694                 /* Check if the current end step is closer to the desired
 695                  * block boundary than the next end step.
 696                  */
 697                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 698                 if (eee[nb].nst > 0 &&
 699                     bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
 700                 {
 701                     set_ee_av(&eee[nb]);
 702                 }
 703                 if (f == 0)
 704                 {
 705                     eee[nb].nst = 1;
 706                 }
 707                 else
 708                 {
 709                     eee[nb].nst += edat->step[f] - edat->step[f-1];
 710                 }
 711                 if (ed->bExactStat)
 712                 {
 713                     add_ee_sum(&eee[nb].sum,es->sum,edat->points[f]);
 714                 }
 715                 else
 716                 {
 717                     add_ee_sum(&eee[nb].sum,edat->s[i].ener[f],1);
 718                 }
 719                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 720                 if (edat->step[f]*nb >= bound_nb)
 721                 {
 722                     set_ee_av(&eee[nb]);
 723                 }
 724             }
 725         }
 726
 727         edat->s[i].av = sum/np;
 728         if (ed->bExactStat)
 729         {
 730             edat->s[i].rmsd = sqrt(sum2/np);
 731         }
 732         else
 733         {
 734             edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
 735         }
 736
 737         if (edat->nframes > 1)
 738         {
 739             edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
 740         }
 741         else
 742         {
 743             edat->s[i].slope = 0;
 744         }
 745
 746         nee  = 0;
 747         see2 = 0;
 748         for(nb=nbmin; nb<=nbmax; nb++)
 749         {
 750             /* Check if we actually got nb blocks and if the smallest
 751              * block is not shorter than 80% of the average.
 752              */
 753             if (debug)
 754             {
 755                 char buf1[STEPSTRSIZE],buf2[STEPSTRSIZE];
 756                 fprintf(debug,"Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
 757                         nb,eee[nb].b,
 758                         gmx_step_str(eee[nb].nst_min,buf1),
 759                         gmx_step_str(edat->nsteps,buf2));
 760             }
 761             if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
 762             {
 763                 see2 += calc_ee2(nb,&eee[nb].sum);
 764                 nee++;
 765             }
 766         }
 767         if (nee > 0)
 768         {
 769             edat->s[i].ee = sqrt(see2/nee);
 770         }
 771         else
 772         {
 773             edat->s[i].ee = -1;
 774         }
 775     }
 776     sfree(eee);
 777 }
 778
 779 static enerdata_t *calc_sum(int nset,enerdata_t *edat,int nbmin,int nbmax)
 780 {
 781     enerdata_t *esum;
 782     enerdat_t *s;
 783     int f,i;
 784     double sum;
 785
 786     snew(esum,1);
 787     *esum = *edat;
 788     snew(esum->s,1);
 789     s = &esum->s[0];
 790     snew(s->ener,esum->nframes);
 791     snew(s->es  ,esum->nframes);
 792
 793     s->bExactStat = TRUE;
 794     s->slope      = 0;
 795     for(i=0; i<nset; i++)
 796     {
 797         if (!edat->s[i].bExactStat)
 798         {
 799             s->bExactStat = FALSE;
 800         }
 801         s->slope += edat->s[i].slope;
 802     }
 803
 804     for(f=0; f<edat->nframes; f++)
 805     {
 806         sum = 0;
 807         for(i=0; i<nset; i++)
 808         {
 809             sum += edat->s[i].ener[f];
 810         }
 811         s->ener[f] = sum;
 812         sum = 0;
 813         for(i=0; i<nset; i++)
 814         {
 815             sum += edat->s[i].es[f].sum;
 816         }
 817         s->es[f].sum  = sum;
 818         s->es[f].sum2 = 0;
 819     }
 820
 821     calc_averages(1,esum,nbmin,nbmax);
 822
 823     return esum;
 824 }
 825
 826 static char *ee_pr(double ee,char *buf)
 827 {
 828     char   tmp[100];
 829     double rnd;
 830
 831     if (ee < 0)
 832     {
 833         sprintf(buf,"%s","--");
 834     }
 835     else
 836     {
 837         /* Round to two decimals by printing. */
 838         sprintf(tmp,"%.1e",ee);
 839         sscanf(tmp,"%lf",&rnd);
 840         sprintf(buf,"%g",rnd);
 841     }
 842
 843     return buf;
 844 }
 845
 846 static void remove_drift(int nset,int nbmin,int nbmax,real dt,enerdata_t *edat)
 847 {
 848 /* Remove the drift by performing a fit to y = ax+b.
 849    Uses 5 iterations. */
 850     int i,j,k;
 851     double delta,d,sd,sd2;
 852
 853     edat->npoints = edat->nframes;
 854     edat->nsteps = edat->nframes;
 855
 856     for(k=0; (k<5); k++)
 857     {
 858         for(i=0; (i<nset); i++)
 859         {
 860             delta = edat->s[i].slope*dt;
 861
 862             if (NULL != debug)
 863                 fprintf(debug,"slope for set %d is %g\n",i,edat->s[i].slope);
 864
 865             for(j=0; (j<edat->nframes); j++)
 866             {
 867                 edat->s[i].ener[j]   -= j*delta;
 868                 edat->s[i].es[j].sum  = 0;
 869                 edat->s[i].es[j].sum2 = 0;
 870             }
 871         }
 872         calc_averages(nset,edat,nbmin,nbmax);
 873     }
 874 }
 875
 876 static void calc_fluctuation_props(FILE *fp,
 877                                    gmx_bool bDriftCorr,real dt,
 878                                    int nset,int set[],int nmol,
 879                                    char **leg,enerdata_t *edat,
 880                                    int nbmin,int nbmax)
 881 {
 882     int i,j;
 883     double vvhh,vv,v,h,hh2,vv2,varv,hh,varh,tt,cv,cp,alpha,kappa,dcp,et,varet;
 884     double NANO3;
 885     enum { eVol, eEnth, eTemp, eEtot, eNR };
 886     const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
 887     int ii[eNR];
 888
 889     NANO3 = NANO*NANO*NANO;
 890     if (!bDriftCorr)
 891     {
 892         fprintf(fp,"\nYou may want to use the -driftcorr flag in order to correct\n"
 893                 "for spurious drift in the graphs. Note that this is not\n"
 894                 "a substitute for proper equilibration and sampling!\n");
 895     }
 896     else
 897     {
 898         remove_drift(nset,nbmin,nbmax,dt,edat);
 899     }
 900     for(i=0; (i<eNR); i++)
 901     {
 902         for(ii[i]=0; (ii[i]<nset &&
 903                       (gmx_strcasecmp(leg[ii[i]],my_ener[i]) != 0)); ii[i]++)
 904             ;
 905 /*        if (ii[i] < nset)
 906             fprintf(fp,"Found %s data.\n",my_ener[i]);
 907 */  }
 908     /* Compute it all! */
 909     vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
 910
 911     /* Temperature */
 912     tt = NOTSET;
 913     if (ii[eTemp] < nset)
 914     {
 915         tt    = edat->s[ii[eTemp]].av;
 916     }
 917     /* Volume */
 918     vv = varv = NOTSET;
 919     if ((ii[eVol] < nset) && (ii[eTemp] < nset))
 920     {
 921         vv    = edat->s[ii[eVol]].av*NANO3;
 922         varv  = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
 923         kappa = (varv/vv)/(BOLTZMANN*tt);
 924     }
 925     /* Enthalpy */
 926     hh = varh = NOTSET;
 927     if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
 928     {
 929         hh    = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
 930         varh  = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
 931         cp    = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
 932     }
 933     /* Total energy */
 934     et = varet = NOTSET;
 935     if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
 936     {
 937         /* Only compute cv in constant volume runs, which we can test
 938            by checking whether the enthalpy was computed.
 939         */
 940         et    = edat->s[ii[eEtot]].av;
 941         varet = sqr(edat->s[ii[eEtot]].rmsd);
 942         cv    = KILO*((varet/nmol)/(BOLTZ*tt*tt));
 943     }
 944     /* Alpha, dcp */
 945     if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
 946     {
 947         vvhh = 0;
 948         for(j=0; (j<edat->nframes); j++)
 949         {
 950             v = edat->s[ii[eVol]].ener[j]*NANO3;
 951             h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
 952             vvhh += (v*h);
 953         }
 954         vvhh /= edat->nframes;
 955         alpha = (vvhh-vv*hh)/(vv*BOLTZMANN*tt*tt);
 956         dcp   = (vv*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
 957     }
 958
 959     if (tt != NOTSET)
 960     {
 961         if (nmol < 2)
 962             fprintf(fp,"\nWARNING: nmol = %d, this may not be what you want.\n",
 963                     nmol);
 964         fprintf(fp,"\nTemperature dependent fluctuation properties at T = %g.\n",tt);
 965         fprintf(fp,"\nHeat capacities obtained from fluctuations do *not* include\n");
 966         fprintf(fp,"quantum corrections. If you want to get a more accurate estimate\n");
 967         fprintf(fp,"please use the g_dos program.\n\n");
 968         fprintf(fp,"WARNING: Please verify that your simulations are converged and perform\n"
 969                 "a block-averaging error analysis (not implemented in g_energy yet)\n");
 970
 971         if (debug != NULL)
 972         {
 973             if (varv != NOTSET)
 974                 fprintf(fp,"varv  =  %10g (m^6)\n",varv*AVOGADRO/nmol);
 975             if (vvhh != NOTSET)
 976                 fprintf(fp,"vvhh  =  %10g (m^3 J)\n",vvhh);
 977         }
 978         if (vv != NOTSET)
 979             fprintf(fp,"Volume                                   = %10g m^3/mol\n",
 980                     vv*AVOGADRO/nmol);
 981         if (varh != NOTSET)
 982             fprintf(fp,"Enthalpy                                 = %10g kJ/mol\n",
 983                 hh*AVOGADRO/(KILO*nmol));
 984         if (alpha != NOTSET)
 985             fprintf(fp,"Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
 986                     alpha);
 987         if (kappa != NOTSET)
 988         {
 989             fprintf(fp,"Isothermal Compressibility Kappa         = %10g (J/m^3)\n",
 990                     kappa);
 991             fprintf(fp,"Adiabatic bulk modulus                   = %10g (m^3/J)\n",
 992                     1.0/kappa);
 993         }
 994         if (cp != NOTSET)
 995             fprintf(fp,"Heat capacity at constant pressure Cp    = %10g J/mol K\n",
 996                     cp);
 997         if (cv != NOTSET)
 998             fprintf(fp,"Heat capacity at constant volume Cv      = %10g J/mol K\n",
 999                     cv);
1000         if (dcp != NOTSET)
1001             fprintf(fp,"Cp-Cv                                    =  %10g J/mol K\n",
1002                     dcp);
1003         please_cite(fp,"Allen1987a");
1004     }
1005     else
1006     {
1007         fprintf(fp,"You should select the temperature in order to obtain fluctuation properties.\n");
1008     }
1009 }
1010
1011 static void analyse_ener(gmx_bool bCorr,const char *corrfn,
1012                          gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,gmx_bool bTempFluct,
1013                          gmx_bool bVisco,const char *visfn,int nmol,
1014                          gmx_large_int_t start_step,double start_t,
1015                          gmx_large_int_t step,double t,
1016                          double time[], real reftemp,
1017                          enerdata_t *edat,
1018                          int nset,int set[],gmx_bool *bIsEner,
1019                          char **leg,gmx_enxnm_t *enm,
1020                          real Vaver,real ezero,
1021                          int nbmin,int nbmax,
1022                          const output_env_t oenv)
1023 {
1024   FILE *fp;
1025   /* Check out the printed manual for equations! */
1026   double Dt,aver,stddev,errest,delta_t,totaldrift;
1027   enerdata_t *esum=NULL;
1028   real xxx,integral,intBulk,Temp=0,Pres=0;
1029   real sfrac,oldfrac,diffsum,diffav,fstep,pr_aver,pr_stddev,pr_errest;
1030   double beta=0,expE,expEtot,*fee=NULL;
1031   gmx_large_int_t nsteps;
1032   int  nexact,nnotexact;
1033   double x1m,x1mk;
1034   int  i,j,k,nout;
1035   real chi2;
1036   char buf[256],eebuf[100];
1037
1038   nsteps  = step - start_step + 1;
1039   if (nsteps < 1) {
1040     fprintf(stdout,"Not enough steps (%s) for statistics\n",
1041             gmx_step_str(nsteps,buf));
1042   }
1043   else {
1044     /* Calculate the time difference */
1045     delta_t = t - start_t;
1046
1047     fprintf(stdout,"\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1048             gmx_step_str(nsteps,buf),start_t,t,nset);
1049
1050     calc_averages(nset,edat,nbmin,nbmax);
1051
1052     if (bSum) {
1053         esum = calc_sum(nset,edat,nbmin,nbmax);
1054     }
1055
1056     if (edat->npoints == 0) {
1057       nexact    = 0;
1058       nnotexact = nset;
1059     } else {
1060       nexact    = 0;
1061       nnotexact = 0;
1062       for(i=0; (i<nset); i++) {
1063         if (edat->s[i].bExactStat) {
1064           nexact++;
1065         } else {
1066           nnotexact++;
1067         }
1068       }
1069     }
1070
1071     if (nnotexact == 0) {
1072       fprintf(stdout,"All statistics are over %s points\n",
1073               gmx_step_str(edat->npoints,buf));
1074     } else if (nexact == 0 || edat->npoints == edat->nframes) {
1075       fprintf(stdout,"All statistics are over %d points (frames)\n",
1076               edat->nframes);
1077     } else {
1078       fprintf(stdout,"The term%s",nnotexact==1 ? "" : "s");
1079       for(i=0; (i<nset); i++) {
1080         if (!edat->s[i].bExactStat) {
1081           fprintf(stdout," '%s'",leg[i]);
1082         }
1083       }
1084       fprintf(stdout," %s has statistics over %d points (frames)\n",
1085               nnotexact==1 ? "is" : "are",edat->nframes);
1086       fprintf(stdout,"All other statistics are over %s points\n",
1087               gmx_step_str(edat->npoints,buf));
1088     }
1089     fprintf(stdout,"\n");
1090
1091     fprintf(stdout,"%-24s %10s %10s %10s %10s",
1092             "Energy","Average","Err.Est.","RMSD","Tot-Drift");
1093     if (bFee)
1094       fprintf(stdout,"  %10s\n","-kT ln<e^(E/kT)>");
1095     else
1096       fprintf(stdout,"\n");
1097     fprintf(stdout,"-------------------------------------------------------------------------------\n");
1098
1099     /* Initiate locals, only used with -sum */
1100     expEtot=0;
1101     if (bFee) {
1102       beta = 1.0/(BOLTZ*reftemp);
1103       snew(fee,nset);
1104     }
1105     for(i=0; (i<nset); i++) {
1106       aver   = edat->s[i].av;
1107       stddev = edat->s[i].rmsd;
1108       errest = edat->s[i].ee;
1109
1110       if (bFee) {
1111         expE = 0;
1112         for(j=0; (j<edat->nframes); j++) {
1113           expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1114         }
1115         if (bSum)
1116           expEtot+=expE/edat->nframes;
1117
1118         fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1119       }
1120       if (strstr(leg[i],"empera") != NULL) {
1121         Temp = aver;
1122       } else if (strstr(leg[i],"olum") != NULL) {
1123         Vaver= aver;
1124       } else if (strstr(leg[i],"essure") != NULL) {
1125         Pres = aver;
1126       }
1127       if (bIsEner[i]) {
1128         pr_aver   = aver/nmol-ezero;
1129         pr_stddev = stddev/nmol;
1130         pr_errest = errest/nmol;
1131       }
1132       else {
1133         pr_aver   = aver;
1134         pr_stddev = stddev;
1135         pr_errest = errest;
1136       }
1137
1138       /* Multiply the slope in steps with the number of steps taken */
1139       totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1140       if (bIsEner[i])
1141       {
1142           totaldrift /= nmol;
1143       }
1144
1145       fprintf(stdout,"%-24s %10g %10s %10g %10g",
1146               leg[i],pr_aver,ee_pr(pr_errest,eebuf),pr_stddev,totaldrift);
1147       if (bFee)
1148         fprintf(stdout,"  %10g",fee[i]);
1149
1150       fprintf(stdout,"  (%s)\n",enm[set[i]].unit);
1151
1152       if (bFluct) {
1153         for(j=0; (j<edat->nframes); j++)
1154           edat->s[i].ener[j] -= aver;
1155       }
1156     }
1157     if (bSum) {
1158         totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1159       fprintf(stdout,"%-24s %10g %10s %10s %10g  (%s)",
1160               "Total",esum->s[0].av/nmol,ee_pr(esum->s[0].ee/nmol,eebuf),
1161               "--",totaldrift/nmol,enm[set[0]].unit);
1162       /* pr_aver,pr_stddev,a,totaldrift */
1163       if (bFee)
1164         fprintf(stdout,"  %10g  %10g\n",
1165                 log(expEtot)/beta + esum->s[0].av/nmol,log(expEtot)/beta);
1166       else
1167         fprintf(stdout,"\n");
1168     }
1169
1170     /* Do correlation function */
1171     if (edat->nframes > 1)
1172     {
1173         Dt = delta_t/(edat->nframes - 1);
1174     }
1175     else
1176     {
1177         Dt = 0;
1178     }
1179     if (bVisco) {
1180       const char* leg[] = { "Shear", "Bulk" };
1181       real factor;
1182       real **eneset;
1183       real **enesum;
1184
1185       /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1186
1187       /* Symmetrise tensor! (and store in first three elements)
1188        * And subtract average pressure!
1189        */
1190       snew(eneset,12);
1191       for(i=0; i<12; i++) {
1192           snew(eneset[i],edat->nframes);
1193       }
1194       snew(enesum,3);
1195       for(i=0; i<3; i++) {
1196         snew(enesum[i],edat->nframes);
1197       }
1198       for(i=0; (i<edat->nframes); i++) {
1199         eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1200         eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1201         eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1202         for(j=3; j<=11; j++) {
1203           eneset[j][i] = edat->s[j].ener[i];
1204         }
1205         eneset[11][i] -= Pres;
1206         enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1207         enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1208         enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1209       }
1210
1211       einstein_visco("evisco.xvg","eviscoi.xvg",
1212                      3,edat->nframes,enesum,Vaver,Temp,nsteps,time,oenv);
1213
1214       /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1215       /* Do it for shear viscosity */
1216       strcpy(buf,"Shear Viscosity");
1217       low_do_autocorr(corrfn,oenv,buf,edat->nframes,3,
1218                       (edat->nframes+1)/2,eneset,Dt,
1219                       eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1220
1221       /* Now for bulk viscosity */
1222       strcpy(buf,"Bulk Viscosity");
1223       low_do_autocorr(corrfn,oenv,buf,edat->nframes,1,
1224                       (edat->nframes+1)/2,&(eneset[11]),Dt,
1225                       eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1226
1227       factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1228       fp=xvgropen(visfn,buf,"Time (ps)","\\8h\\4 (cp)",oenv);
1229       xvgr_legend(fp,asize(leg),leg,oenv);
1230
1231       /* Use trapezium rule for integration */
1232       integral = 0;
1233       intBulk  = 0;
1234       nout = get_acfnout();
1235       if ((nout < 2) || (nout >= edat->nframes/2))
1236           nout = edat->nframes/2;
1237       for(i=1; (i<nout); i++)
1238       {
1239           integral += 0.5*(eneset[0][i-1]  + eneset[0][i])*factor;
1240           intBulk  += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1241           fprintf(fp,"%10g  %10g  %10g\n",(i*Dt),integral,intBulk);
1242       }
1243       ffclose(fp);
1244     }
1245     else if (bCorr) {
1246       if (bFluct)
1247           strcpy(buf,"Autocorrelation of Energy Fluctuations");
1248       else
1249           strcpy(buf,"Energy Autocorrelation");
1250 #if 0
1251       do_autocorr(corrfn,oenv,buf,edat->nframes,
1252                   bSum ? 1                 : nset,
1253                   bSum ? &edat->s[nset-1].ener : eneset,
1254                   (delta_t/edat->nframes),eacNormal,FALSE);
1255 #endif
1256     }
1257   }
1258 }
1259
1260 static void print_time(FILE *fp,double t)
1261 {
1262   fprintf(fp,"%12.6f",t);
1263 }
1264
1265 static void print1(FILE *fp,gmx_bool bDp,real e)
1266 {
1267   if (bDp)
1268     fprintf(fp,"  %16.12f",e);
1269   else
1270     fprintf(fp,"  %10.6f",e);
1271 }
1272
1273 static void fec(const char *ene2fn, const char *runavgfn,
1274                 real reftemp, int nset, int set[], char *leg[],
1275                 enerdata_t *edat, double time[],
1276                 const output_env_t oenv)
1277 {
1278   const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1279                            "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1280   FILE *fp;
1281   ener_file_t enx;
1282   int  nre,timecheck,step,nenergy,nenergy2,maxenergy;
1283   int  i,j;
1284   gmx_bool bCont;
1285   real aver, beta;
1286   real **eneset2;
1287   double dE, sum;
1288   gmx_enxnm_t *enm=NULL;
1289   t_enxframe *fr;
1290   char buf[22];
1291
1292   /* read second energy file */
1293   snew(fr,1);
1294   enm = NULL;
1295   enx = open_enx(ene2fn,"r");
1296   do_enxnms(enx,&(fr->nre),&enm);
1297
1298   snew(eneset2,nset+1);
1299   nenergy2=0;
1300   maxenergy=0;
1301   timecheck=0;
1302   do {
1303     /* This loop searches for the first frame (when -b option is given),
1304      * or when this has been found it reads just one energy frame
1305      */
1306     do {
1307       bCont = do_enx(enx,fr);
1308
1309       if (bCont)
1310         timecheck = check_times(fr->t);
1311
1312     } while (bCont && (timecheck < 0));
1313
1314     /* Store energies for analysis afterwards... */
1315     if ((timecheck == 0) && bCont) {
1316       if (fr->nre > 0) {
1317         if ( nenergy2 >= maxenergy ) {
1318           maxenergy += 1000;
1319           for(i=0; i<=nset; i++)
1320             srenew(eneset2[i],maxenergy);
1321         }
1322         if (fr->t != time[nenergy2])
1323           fprintf(stderr,"\nWARNING time mismatch %g!=%g at frame %s\n",
1324                   fr->t, time[nenergy2], gmx_step_str(fr->step,buf));
1325         for(i=0; i<nset; i++)
1326           eneset2[i][nenergy2] = fr->ener[set[i]].e;
1327         nenergy2++;
1328       }
1329     }
1330   } while (bCont && (timecheck == 0));
1331
1332   /* check */
1333   if (edat->nframes != nenergy2) {
1334     fprintf(stderr,"\nWARNING file length mismatch %d!=%d\n",
1335             edat->nframes,nenergy2);
1336   }
1337   nenergy = min(edat->nframes,nenergy2);
1338
1339   /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1340   fp=NULL;
1341   if (runavgfn) {
1342     fp=xvgropen(runavgfn,"Running average free energy difference",
1343                 "Time (" unit_time ")","\\8D\\4E (" unit_energy ")",oenv);
1344     xvgr_legend(fp,asize(ravgleg),ravgleg,oenv);
1345   }
1346   fprintf(stdout,"\n%-24s %10s\n",
1347           "Energy","dF = -kT ln < exp(-(EB-EA)/kT) >A");
1348   sum=0;
1349   beta = 1.0/(BOLTZ*reftemp);
1350   for(i=0; i<nset; i++) {
1351     if (gmx_strcasecmp(leg[i],enm[set[i]].name)!=0)
1352       fprintf(stderr,"\nWARNING energy set name mismatch %s!=%s\n",
1353               leg[i],enm[set[i]].name);
1354     for(j=0; j<nenergy; j++) {
1355       dE = eneset2[i][j] - edat->s[i].ener[j];
1356       sum += exp(-dE*beta);
1357       if (fp)
1358         fprintf(fp,"%10g %10g %10g\n",
1359                 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1360     }
1361     aver = -BOLTZ*reftemp*log(sum/nenergy);
1362     fprintf(stdout,"%-24s %10g\n",leg[i],aver);
1363   }
1364   if(fp) ffclose(fp);
1365   sfree(fr);
1366 }
1367
1368
1369 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, const char *filename, gmx_bool bDp,
1370                     int *blocks, int *hists, int *samples, int *nlambdas, const output_env_t oenv)
1371 {
1372     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
1373     char title[STRLEN],label_x[STRLEN],label_y[STRLEN], legend[STRLEN];
1374     char buf[STRLEN];
1375     gmx_bool first=FALSE;
1376     int nblock_hist=0, nblock_dh=0, nblock_dhcoll=0;
1377     int i,j,k;
1378     /* coll data */
1379     double temp=0, start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
1380     static int setnr=0;
1381     gmx_bool changing_lambda=FALSE;
1382
1383     /* now count the blocks & handle the global dh data */
1384     for(i=0;i<fr->nblock;i++)
1385     {
1386         if (fr->block[i].id == enxDHHIST)
1387         {
1388             nblock_hist++;
1389         }
1390         else if (fr->block[i].id == enxDH)
1391         {
1392             nblock_dh++;
1393         }
1394         else if (fr->block[i].id == enxDHCOLL)
1395         {
1396             nblock_dhcoll++;
1397             if ( (fr->block[i].nsub < 1) ||
1398                  (fr->block[i].sub[0].type != xdr_datatype_double) ||
1399                  (fr->block[i].sub[0].nr < 5))
1400             {
1401                 gmx_fatal(FARGS, "Unexpected block data");
1402             }
1403
1404             /* read the data from the DHCOLL block */
1405             temp =         fr->block[i].sub[0].dval[0];
1406             start_time =   fr->block[i].sub[0].dval[1];
1407             delta_time =   fr->block[i].sub[0].dval[2];
1408             start_lambda = fr->block[i].sub[0].dval[3];
1409             delta_lambda = fr->block[i].sub[0].dval[4];
1410             changing_lambda = (delta_lambda != 0);
1411         }
1412     }
1413
1414     if (nblock_hist == 0 && nblock_dh == 0)
1415     {
1416         /* don't do anything */
1417         return;
1418     }
1419     if (nblock_hist>0 && nblock_dh>0)
1420     {
1421         gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1422     }
1423     if (! *fp_dhdl )
1424     {
1425         if (nblock_dh>0)
1426         {
1427             /* we have standard, non-histogram data -- call open_dhdl to open the file */
1428             *fp_dhdl=open_dhdl(filename,ir,oenv);
1429         }
1430         else
1431         {
1432             sprintf(title,"N(%s)",deltag);
1433             sprintf(label_x,"%s (%s)",deltag,unit_energy);
1434             sprintf(label_y,"Samples");
1435             *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,oenv);
1436             sprintf(buf,"T = %g (K), %s = %g", temp, lambda, start_lambda);
1437             xvgr_subtitle(*fp_dhdl,buf,oenv);
1438         }
1439     }
1440
1441     (*hists)+=nblock_hist;
1442     (*blocks)+=nblock_dh;
1443     (*nlambdas) = nblock_hist+nblock_dh;
1444
1445     /* write the data */
1446     if (nblock_hist > 0)
1447     {
1448         gmx_large_int_t sum=0;
1449         /* histograms */
1450         for(i=0;i<fr->nblock;i++)
1451         {
1452             t_enxblock *blk=&(fr->block[i]);
1453             if (blk->id==enxDHHIST)
1454             {
1455                 double foreign_lambda, dx;
1456                 gmx_large_int_t x0;
1457                 int nhist, derivative;
1458
1459                 /* check the block types etc. */
1460                 if ( (blk->nsub < 2) ||
1461                      (blk->sub[0].type != xdr_datatype_double) ||
1462                      (blk->sub[1].type != xdr_datatype_large_int) ||
1463                      (blk->sub[0].nr < 2)  ||
1464                      (blk->sub[1].nr < 2) )
1465                 {
1466                     gmx_fatal(FARGS, "Unexpected block data in file");
1467                 }
1468                 foreign_lambda=blk->sub[0].dval[0];
1469                 dx=blk->sub[0].dval[1];
1470                 nhist=blk->sub[1].lval[0];
1471                 derivative=blk->sub[1].lval[1];
1472                 for(j=0;j<nhist;j++)
1473                 {
1474                     const char *lg[1];
1475                     x0=blk->sub[1].lval[2+j];
1476
1477                     if (!derivative)
1478                     {
1479                         sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1480                                 deltag, lambda, foreign_lambda,
1481                                 lambda, start_lambda);
1482                     }
1483                     else
1484                     {
1485                         sprintf(legend, "N(%s | %s=%g)",
1486                                 dhdl, lambda, start_lambda);
1487                     }
1488
1489                     lg[0]=legend;
1490                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1491                     setnr++;
1492                     for(k=0;k<blk->sub[j+2].nr;k++)
1493                     {
1494                         int hist;
1495                         double xmin, xmax;
1496
1497                         hist=blk->sub[j+2].ival[k];
1498                         xmin=(x0+k)*dx;
1499                         xmax=(x0+k+1)*dx;
1500                         fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
1501                                 xmax, hist);
1502                         sum+=hist;
1503                     }
1504                     /* multiple histogram data blocks in one histogram
1505                        mean that the second one is the reverse of the first one:
1506                        for dhdl derivatives, it's important to know both the
1507                        maximum and minimum values */
1508                     dx=-dx;
1509                 }
1510             }
1511         }
1512         (*samples) += (int)(sum/nblock_hist);
1513     }
1514     else
1515     {
1516         /* raw dh */
1517         int len=0;
1518         char **setnames=NULL;
1519         int nnames=nblock_dh;
1520
1521         for(i=0;i<fr->nblock;i++)
1522         {
1523             t_enxblock *blk=&(fr->block[i]);
1524             if (blk->id == enxDH)
1525             {
1526                 if (len == 0)
1527                 {
1528                     len=blk->sub[2].nr;
1529                 }
1530                 else
1531                 {
1532                     if (len!=blk->sub[2].nr)
1533                     {
1534                         gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1535                     }
1536                 }
1537             }
1538         }
1539         (*samples) += len;
1540
1541         for(i=0;i<len;i++)
1542         {
1543             double time=start_time + delta_time*i;
1544
1545             fprintf(*fp_dhdl,"%.4f ", time);
1546
1547             for(j=0;j<fr->nblock;j++)
1548             {
1549                 t_enxblock *blk=&(fr->block[j]);
1550                 if (blk->id == enxDH)
1551                 {
1552                     double value;
1553                     if (blk->sub[2].type == xdr_datatype_float)
1554                     {
1555                         value=blk->sub[2].fval[i];
1556                     }
1557                     else
1558                     {
1559                         value=blk->sub[2].dval[i];
1560                     }
1561                     /* we need to decide which data type it is based on the count*/
1562
1563                     if (j==1 && ir->bExpanded)
1564                     {
1565                         fprintf(*fp_dhdl,"%4d", (int)value);   /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
1566                     } else {
1567                         if (bDp) {
1568                             fprintf(*fp_dhdl," %#.12g", value);   /* print normal precision */
1569                         }
1570                         else
1571                         {
1572                             fprintf(*fp_dhdl," %#.8g", value);   /* print normal precision */
1573                         }
1574                     }
1575                 }
1576             }
1577             fprintf(*fp_dhdl, "\n");
1578         }
1579     }
1580 }
1581
1582
1583 int gmx_energy(int argc,char *argv[])
1584 {
1585   const char *desc[] = {
1586
1587     "[TT]g_energy[tt] extracts energy components or distance restraint",
1588     "data from an energy file. The user is prompted to interactively",
1589     "select the desired energy terms.[PAR]",
1590
1591     "Average, RMSD, and drift are calculated with full precision from the",
1592     "simulation (see printed manual). Drift is calculated by performing",
1593     "a least-squares fit of the data to a straight line. The reported total drift",
1594     "is the difference of the fit at the first and last point.",
1595     "An error estimate of the average is given based on a block averages",
1596     "over 5 blocks using the full-precision averages. The error estimate",
1597     "can be performed over multiple block lengths with the options",
1598     "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1599     "[BB]Note[bb] that in most cases the energy files contains averages over all",
1600     "MD steps, or over many more points than the number of frames in",
1601     "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
1602     "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1603     "file, the statistics mentioned above are simply over the single, per-frame",
1604     "energy values.[PAR]",
1605
1606     "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1607
1608     "Some fluctuation-dependent properties can be calculated provided",
1609     "the correct energy terms are selected. The following properties",
1610     "will be computed:[BR]",
1611     "Property                        Energy terms needed[BR]",
1612     "---------------------------------------------------[BR]",
1613     "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp     [BR]",
1614     "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp         [BR]",
1615     "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1616     "Isothermal compressibility:     Vol, Temp          [BR]",
1617     "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1618     "---------------------------------------------------[BR]",
1619     "You always need to set the number of molecules [TT]-nmol[tt].",
1620     "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1621     "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
1622     "When the [TT]-viol[tt] option is set, the time averaged",
1623     "violations are plotted and the running time-averaged and",
1624     "instantaneous sum of violations are recalculated. Additionally",
1625     "running time-averaged and instantaneous distances between",
1626     "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1627
1628     "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1629     "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1630     "The first two options plot the orientation, the last three the",
1631     "deviations of the orientations from the experimental values.",
1632     "The options that end on an 'a' plot the average over time",
1633     "as a function of restraint. The options that end on a 't'",
1634     "prompt the user for restraint label numbers and plot the data",
1635     "as a function of time. Option [TT]-odr[tt] plots the RMS",
1636     "deviation as a function of restraint.",
1637     "When the run used time or ensemble averaged orientation restraints,",
1638     "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1639     "not ensemble-averaged orientations and deviations instead of",
1640     "the time and ensemble averages.[PAR]",
1641
1642     "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1643     "tensor for each orientation restraint experiment. With option",
1644     "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1645
1646     "Option [TT]-odh[tt] extracts and plots the free energy data",
1647     "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1648     "from the [TT]ener.edr[tt] file.[PAR]",
1649
1650     "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1651     "difference with an ideal gas state: [BR]",
1652     "  [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
1653     "  [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
1654     "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1655     "the average is over the ensemble (or time in a trajectory).",
1656     "Note that this is in principle",
1657     "only correct when averaging over the whole (Boltzmann) ensemble",
1658     "and using the potential energy. This also allows for an entropy",
1659     "estimate using:[BR]",
1660     "  [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]",
1661     "  [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
1662     "[PAR]",
1663
1664     "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1665     "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1666     "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1667     "files, and the average is over the ensemble A. The running average",
1668     "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1669     "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1670
1671   };
1672   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
1673   static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
1674   static int  skip=0,nmol=1,nbmin=5,nbmax=5;
1675   static real reftemp=300.0,ezero=0;
1676   t_pargs pa[] = {
1677     { "-fee",   FALSE, etBOOL,  {&bFee},
1678       "Do a free energy estimate" },
1679     { "-fetemp", FALSE, etREAL,{&reftemp},
1680       "Reference temperature for free energy calculation" },
1681     { "-zero", FALSE, etREAL, {&ezero},
1682       "Subtract a zero-point energy" },
1683     { "-sum",  FALSE, etBOOL, {&bSum},
1684       "Sum the energy terms selected rather than display them all" },
1685     { "-dp",   FALSE, etBOOL, {&bDp},
1686       "Print energies in high precision" },
1687     { "-nbmin", FALSE, etINT, {&nbmin},
1688       "Minimum number of blocks for error estimate" },
1689     { "-nbmax", FALSE, etINT, {&nbmax},
1690       "Maximum number of blocks for error estimate" },
1691     { "-mutot",FALSE, etBOOL, {&bMutot},
1692       "Compute the total dipole moment from the components" },
1693     { "-skip", FALSE, etINT,  {&skip},
1694       "Skip number of frames between data points" },
1695     { "-aver", FALSE, etBOOL, {&bPrAll},
1696       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1697     { "-nmol", FALSE, etINT,  {&nmol},
1698       "Number of molecules in your sample: the energies are divided by this number" },
1699     { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1700       "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1701     { "-fluc", FALSE, etBOOL, {&bFluct},
1702       "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1703     { "-orinst", FALSE, etBOOL, {&bOrinst},
1704       "Analyse instantaneous orientation data" },
1705     { "-ovec", FALSE, etBOOL, {&bOvec},
1706       "Also plot the eigenvectors with [TT]-oten[tt]" }
1707   };
1708   const char* drleg[] = {
1709     "Running average",
1710     "Instantaneous"
1711   };
1712   static const char *setnm[] = {
1713     "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1714     "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1715     "Volume",  "Pressure"
1716   };
1717
1718   FILE       *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1719   FILE       *fp_dhdl=NULL;
1720   FILE       **drout;
1721   ener_file_t fp;
1722   int        timecheck=0;
1723   gmx_mtop_t mtop;
1724   gmx_localtop_t *top=NULL;
1725   t_inputrec ir;
1726   t_energy   **ee;
1727   enerdata_t edat;
1728   gmx_enxnm_t *enm=NULL;
1729   t_enxframe *frame,*fr=NULL;
1730   int        cur=0;
1731 #define NEXT (1-cur)
1732   int        nre,teller,teller_disre,nfr;
1733   gmx_large_int_t start_step;
1734   int        nor=0,nex=0,norfr=0,enx_i=0;
1735   real       start_t;
1736   real       *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1737   int        *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1738   int        nbounds=0,npairs;
1739   gmx_bool       bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN,bDHDL;
1740   gmx_bool       bFoundStart,bCont,bEDR,bVisco;
1741   double     sum,sumaver,sumt,ener,dbl;
1742   double     *time=NULL;
1743   real       Vaver;
1744   int        *set=NULL,i,j,k,nset,sss;
1745   gmx_bool       *bIsEner=NULL;
1746   char       **pairleg,**odtleg,**otenleg;
1747   char       **leg=NULL;
1748   char       **nms;
1749   char       *anm_j,*anm_k,*resnm_j,*resnm_k;
1750   int        resnr_j,resnr_k;
1751   const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1752   char       buf[256];
1753   output_env_t oenv;
1754   t_enxblock *blk=NULL;
1755   t_enxblock *blk_disre=NULL;
1756   int        ndisre=0;
1757   int        dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas=0;
1758
1759   t_filenm   fnm[] = {
1760     { efEDR, "-f",    NULL,      ffREAD  },
1761     { efEDR, "-f2",   NULL,      ffOPTRD },
1762     { efTPX, "-s",    NULL,      ffOPTRD },
1763     { efXVG, "-o",    "energy",  ffWRITE },
1764     { efXVG, "-viol", "violaver",ffOPTWR },
1765     { efXVG, "-pairs","pairs",   ffOPTWR },
1766     { efXVG, "-ora",  "orienta", ffOPTWR },
1767     { efXVG, "-ort",  "orientt", ffOPTWR },
1768     { efXVG, "-oda",  "orideva", ffOPTWR },
1769     { efXVG, "-odr",  "oridevr", ffOPTWR },
1770     { efXVG, "-odt",  "oridevt", ffOPTWR },
1771     { efXVG, "-oten", "oriten",  ffOPTWR },
1772     { efXVG, "-corr", "enecorr", ffOPTWR },
1773     { efXVG, "-vis",  "visco",   ffOPTWR },
1774     { efXVG, "-ravg", "runavgdf",ffOPTWR },
1775     { efXVG, "-odh",  "dhdl"    ,ffOPTWR }
1776   };
1777 #define NFILE asize(fnm)
1778   int     npargs;
1779   t_pargs *ppa;
1780
1781   CopyRight(stderr,argv[0]);
1782   npargs = asize(pa);
1783   ppa    = add_acf_pargs(&npargs,pa);
1784   parse_common_args(&argc,argv,
1785                     PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1786                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1787
1788   bDRAll = opt2bSet("-pairs",NFILE,fnm);
1789   bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1790   bORA   = opt2bSet("-ora",NFILE,fnm);
1791   bORT   = opt2bSet("-ort",NFILE,fnm);
1792   bODA   = opt2bSet("-oda",NFILE,fnm);
1793   bODR   = opt2bSet("-odr",NFILE,fnm);
1794   bODT   = opt2bSet("-odt",NFILE,fnm);
1795   bORIRE = bORA || bORT || bODA || bODR || bODT;
1796   bOTEN  = opt2bSet("-oten",NFILE,fnm);
1797   bDHDL  = opt2bSet("-odh",NFILE,fnm);
1798
1799   nset = 0;
1800
1801   snew(frame,2);
1802   fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1803   do_enxnms(fp,&nre,&enm);
1804
1805   Vaver = -1;
1806
1807   bVisco = opt2bSet("-vis",NFILE,fnm);
1808
1809   if ((!bDisRe) && (!bDHDL))
1810   {
1811       if (bVisco) {
1812           nset=asize(setnm);
1813           snew(set,nset);
1814           /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1815           for(j=0; j<nset; j++) {
1816               for(i=0; i<nre; i++) {
1817                   if (strstr(enm[i].name,setnm[j])) {
1818                       set[j]=i;
1819                       break;
1820                   }
1821               }
1822               if (i == nre) {
1823                   if (gmx_strcasecmp(setnm[j],"Volume")==0) {
1824                       printf("Enter the box volume (" unit_volume "): ");
1825                       if(1 != scanf("%lf",&dbl))
1826                       {
1827                           gmx_fatal(FARGS,"Error reading user input");
1828                       }
1829                       Vaver = dbl;
1830                   } else
1831                       gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1832                                 setnm[j]);
1833               }
1834           }
1835       }
1836       else
1837       {
1838           set=select_by_name(nre,enm,&nset);
1839       }
1840       /* Print all the different units once */
1841       sprintf(buf,"(%s)",enm[set[0]].unit);
1842       for(i=1; i<nset; i++) {
1843           for(j=0; j<i; j++) {
1844               if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1845                   break;
1846               }
1847           }
1848           if (j == i) {
1849               strcat(buf,", (");
1850               strcat(buf,enm[set[i]].unit);
1851               strcat(buf,")");
1852           }
1853       }
1854       out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1855                    oenv);
1856
1857       snew(leg,nset+1);
1858       for(i=0; (i<nset); i++)
1859           leg[i] = enm[set[i]].name;
1860       if (bSum) {
1861           leg[nset]=strdup("Sum");
1862           xvgr_legend(out,nset+1,(const char**)leg,oenv);
1863       }
1864       else
1865           xvgr_legend(out,nset,(const char**)leg,oenv);
1866
1867       snew(bIsEner,nset);
1868       for(i=0; (i<nset); i++) {
1869           bIsEner[i] = FALSE;
1870           for (j=0; (j <= F_ETOT); j++)
1871               bIsEner[i] = bIsEner[i] ||
1872                         (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1873       }
1874
1875       if (bPrAll && nset > 1) {
1876           gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1877       }
1878
1879       time = NULL;
1880
1881       if (bORIRE || bOTEN)
1882           get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1883
1884       if (bORIRE) {
1885           if (bOrinst)
1886               enx_i = enxORI;
1887           else
1888               enx_i = enxOR;
1889
1890           if (bORA || bODA)
1891               snew(orient,nor);
1892           if (bODR)
1893               snew(odrms,nor);
1894           if (bORT || bODT) {
1895               fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1896               fprintf(stderr,"End your selection with 0\n");
1897               j = -1;
1898               orsel = NULL;
1899               do {
1900                   j++;
1901                   srenew(orsel,j+1);
1902                   if(1 != scanf("%d",&(orsel[j])))
1903                   {
1904                       gmx_fatal(FARGS,"Error reading user input");
1905                   }
1906               } while (orsel[j] > 0);
1907               if (orsel[0] == -1) {
1908                   fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1909                   norsel = nor;
1910                   srenew(orsel,nor);
1911                   for(i=0; i<nor; i++)
1912                       orsel[i] = i;
1913               } else {
1914                   /* Build the selection */
1915                   norsel=0;
1916                   for(i=0; i<j; i++) {
1917                       for(k=0; k<nor; k++)
1918                           if (or_label[k] == orsel[i]) {
1919                               orsel[norsel] = k;
1920                               norsel++;
1921                               break;
1922                           }
1923                       if (k == nor)
1924                           fprintf(stderr,"Orientation restraint label %d not found\n",
1925                                   orsel[i]);
1926                   }
1927               }
1928               snew(odtleg,norsel);
1929               for(i=0; i<norsel; i++) {
1930                   snew(odtleg[i],256);
1931                   sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1932               }
1933               if (bORT) {
1934                   fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1935                                 "Time (ps)","",oenv);
1936                   if (bOrinst)
1937                       fprintf(fort,"%s",orinst_sub);
1938                   xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
1939               }
1940               if (bODT) {
1941                   fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1942                                 "Orientation restraint deviation",
1943                                 "Time (ps)","",oenv);
1944                   if (bOrinst)
1945                       fprintf(fodt,"%s",orinst_sub);
1946                   xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
1947               }
1948           }
1949       }
1950       if (bOTEN) {
1951           foten=xvgropen(opt2fn("-oten",NFILE,fnm),
1952                          "Order tensor","Time (ps)","",oenv);
1953           snew(otenleg,bOvec ? nex*12 : nex*3);
1954           for(i=0; i<nex; i++) {
1955               for(j=0; j<3; j++) {
1956                   sprintf(buf,"eig%d",j+1);
1957                   otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
1958               }
1959               if (bOvec) {
1960                   for(j=0; j<9; j++) {
1961                       sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
1962                       otenleg[12*i+3+j] = strdup(buf);
1963                   }
1964               }
1965           }
1966           xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
1967       }
1968   }
1969   else if (bDisRe)
1970   {
1971       nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
1972                          &mtop,&top,&ir);
1973       snew(violaver,npairs);
1974       out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
1975                    "Time (ps)","nm",oenv);
1976       xvgr_legend(out,2,drleg,oenv);
1977       if (bDRAll) {
1978           fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
1979                             "Time (ps)","Distance (nm)",oenv);
1980           if (output_env_get_print_xvgr_codes(oenv))
1981               fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
1982                       ir.dr_tau);
1983       }
1984   } else if (bDHDL) {
1985       get_dhdl_parms(ftp2fn(efTPX,NFILE,fnm),&ir);
1986   }
1987
1988    /* Initiate energies and set them to zero */
1989    edat.nsteps  = 0;
1990    edat.npoints = 0;
1991    edat.nframes = 0;
1992    edat.step    = NULL;
1993    edat.steps   = NULL;
1994    edat.points  = NULL;
1995    snew(edat.s,nset);
1996
1997    /* Initiate counters */
1998    teller       = 0;
1999    teller_disre = 0;
2000    bFoundStart  = FALSE;
2001    start_step   = 0;
2002    start_t      = 0;
2003    do {
2004      /* This loop searches for the first frame (when -b option is given),
2005       * or when this has been found it reads just one energy frame
2006       */
2007      do {
2008          bCont = do_enx(fp,&(frame[NEXT]));
2009          if (bCont) {
2010              timecheck = check_times(frame[NEXT].t);
2011          }
2012      } while (bCont && (timecheck < 0));
2013
2014      if ((timecheck == 0) && bCont) {
2015        /* We read a valid frame, so we can use it */
2016        fr = &(frame[NEXT]);
2017
2018        if (fr->nre > 0) {
2019          /* The frame contains energies, so update cur */
2020          cur  = NEXT;
2021
2022              if (edat.nframes % 1000 == 0)
2023              {
2024                  srenew(edat.step,edat.nframes+1000);
2025                  memset(&(edat.step[edat.nframes]),0,1000*sizeof(edat.step[0]));
2026                  srenew(edat.steps,edat.nframes+1000);
2027                  memset(&(edat.steps[edat.nframes]),0,1000*sizeof(edat.steps[0]));
2028                  srenew(edat.points,edat.nframes+1000);
2029                  memset(&(edat.points[edat.nframes]),0,1000*sizeof(edat.points[0]));
2030
2031                  for(i=0; i<nset; i++)
2032                  {
2033                      srenew(edat.s[i].ener,edat.nframes+1000);
2034                      memset(&(edat.s[i].ener[edat.nframes]),0,
2035                             1000*sizeof(edat.s[i].ener[0]));
2036                      srenew(edat.s[i].es  ,edat.nframes+1000);
2037                      memset(&(edat.s[i].es[edat.nframes]),0,
2038                             1000*sizeof(edat.s[i].es[0]));
2039                  }
2040              }
2041
2042              nfr = edat.nframes;
2043              edat.step[nfr] = fr->step;
2044
2045              if (!bFoundStart)
2046              {
2047                  bFoundStart = TRUE;
2048                  /* Initiate the previous step data */
2049                  start_step = fr->step;
2050                  start_t    = fr->t;
2051                  /* Initiate the energy sums */
2052                  edat.steps[nfr]  = 1;
2053                  edat.points[nfr] = 1;
2054                  for(i=0; i<nset; i++)
2055                  {
2056                      sss = set[i];
2057                      edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2058                      edat.s[i].es[nfr].sum2 = 0;
2059                  }
2060                  edat.nsteps  = 1;
2061                  edat.npoints = 1;
2062              }
2063              else
2064              {
2065                  edat.steps[nfr] = fr->nsteps;
2066                  {
2067                      if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2068                      {
2069                          if (fr->nsum <= 1)
2070                          {
2071                              edat.points[nfr] = 1;
2072                              for(i=0; i<nset; i++)
2073                              {
2074                                  sss = set[i];
2075                                  edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2076                                  edat.s[i].es[nfr].sum2 = 0;
2077                              }
2078                              edat.npoints += 1;
2079                          }
2080                          else
2081                          {
2082                              edat.points[nfr] = fr->nsum;
2083                              for(i=0; i<nset; i++)
2084                              {
2085                                  sss = set[i];
2086                                  edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2087                                  edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2088                              }
2089                              edat.npoints += fr->nsum;
2090                          }
2091                      }
2092                      else
2093                      {
2094                          /* The interval does not match fr->nsteps:
2095                           * can not do exact averages.
2096                           */
2097                          edat.npoints = 0;
2098                      }
2099                      edat.nsteps = fr->step - start_step + 1;
2100                  }
2101              }
2102              for(i=0; i<nset; i++)
2103              {
2104                  edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2105              }
2106        }
2107        /*
2108         * Define distance restraint legends. Can only be done after
2109         * the first frame has been read... (Then we know how many there are)
2110         */
2111        blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
2112        if (bDisRe && bDRAll && !leg && blk_disre)
2113        {
2114            t_iatom   *fa;
2115            t_iparams *ip;
2116
2117            fa = top->idef.il[F_DISRES].iatoms;
2118            ip = top->idef.iparams;
2119            if (blk_disre->nsub != 2 ||
2120                (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2121            {
2122                gmx_incons("Number of disre sub-blocks not equal to 2");
2123            }
2124
2125            ndisre=blk_disre->sub[0].nr ;
2126            if (ndisre != top->idef.il[F_DISRES].nr/3)
2127            {
2128                gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2129                          ndisre,top->idef.il[F_DISRES].nr/3);
2130            }
2131            snew(pairleg,ndisre);
2132            for(i=0; i<ndisre; i++)
2133            {
2134                snew(pairleg[i],30);
2135                j=fa[3*i+1];
2136                k=fa[3*i+2];
2137                gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
2138                gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
2139                sprintf(pairleg[i],"%d %s %d %s (%d)",
2140                        resnr_j,anm_j,resnr_k,anm_k,
2141                        ip[fa[3*i]].disres.label);
2142            }
2143            set=select_it(ndisre,pairleg,&nset);
2144            snew(leg,2*nset);
2145            for(i=0; (i<nset); i++)
2146            {
2147                snew(leg[2*i],32);
2148                sprintf(leg[2*i],  "a %s",pairleg[set[i]]);
2149                snew(leg[2*i+1],32);
2150                sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
2151            }
2152            xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
2153        }
2154
2155        /*
2156         * Store energies for analysis afterwards...
2157         */
2158        if (!bDisRe && !bDHDL && (fr->nre > 0)) {
2159            if (edat.nframes % 1000 == 0) {
2160                srenew(time,edat.nframes+1000);
2161            }
2162            time[edat.nframes] = fr->t;
2163            edat.nframes++;
2164        }
2165        /*
2166         * Printing time, only when we do not want to skip frames
2167         */
2168        if (!skip || teller % skip == 0) {
2169      if (bDisRe) {
2170        /*******************************************
2171         * D I S T A N C E   R E S T R A I N T S
2172         *******************************************/
2173        if (ndisre > 0)
2174            {
2175  #ifndef GMX_DOUBLE
2176              float *disre_rt =     blk_disre->sub[0].fval;
2177              float *disre_rm3tav = blk_disre->sub[1].fval;
2178  #else
2179              double *disre_rt =     blk_disre->sub[0].dval;
2180              double *disre_rm3tav = blk_disre->sub[1].dval;
2181  #endif
2182
2183          print_time(out,fr->t);
2184          if (violaver == NULL)
2185            snew(violaver,ndisre);
2186
2187          /* Subtract bounds from distances, to calculate violations */
2188          calc_violations(disre_rt, disre_rm3tav,
2189                  nbounds,pair,bounds,violaver,&sumt,&sumaver);
2190
2191          fprintf(out,"  %8.4f  %8.4f\n",sumaver,sumt);
2192          if (bDRAll) {
2193            print_time(fp_pairs,fr->t);
2194            for(i=0; (i<nset); i++) {
2195          sss=set[i];
2196          fprintf(fp_pairs,"  %8.4f", mypow(disre_rm3tav[sss],minthird));
2197          fprintf(fp_pairs,"  %8.4f", disre_rt[sss]);
2198            }
2199            fprintf(fp_pairs,"\n");
2200          }
2201          teller_disre++;
2202        }
2203      }
2204      else if (bDHDL)
2205      {
2206          do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh",NFILE,fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2207      }
2208
2209      /*******************************************
2210       * E N E R G I E S
2211       *******************************************/
2212      else {
2213          if (fr->nre > 0) {
2214              if (bPrAll)
2215              {
2216                  /* We skip frames with single points (usually only the first frame),
2217                   * since they would result in an average plot with outliers.
2218                   */
2219                  if (fr->nsum > 1) {
2220                      print_time(out,fr->t);
2221                       print1(out,bDp,fr->ener[set[0]].e);
2222                       print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
2223                       print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
2224                       fprintf(out,"\n");
2225                  }
2226              }
2227              else
2228              {
2229                  print_time(out,fr->t);
2230                  if (bSum)
2231                  {
2232                      sum = 0;
2233                      for(i=0; i<nset; i++)
2234                      {
2235                          sum += fr->ener[set[i]].e;
2236                      }
2237                      print1(out,bDp,sum/nmol-ezero);
2238                  }
2239                  else
2240                  {
2241                      for(i=0; (i<nset); i++)
2242                      {
2243                          if (bIsEner[i])
2244                          {
2245                              print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
2246                          }
2247                          else
2248                          {
2249                              print1(out,bDp,fr->ener[set[i]].e);
2250                          }
2251                      }
2252                  }
2253                  fprintf(out,"\n");
2254              }
2255          }
2256           blk = find_block_id_enxframe(fr, enx_i, NULL);
2257           if (bORIRE && blk)
2258           {
2259 #ifndef GMX_DOUBLE
2260               xdr_datatype dt=xdr_datatype_float;
2261 #else
2262               xdr_datatype dt=xdr_datatype_double;
2263 #endif
2264               real *vals;
2265
2266               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2267               {
2268                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2269               }
2270 #ifndef GMX_DOUBLE
2271               vals=blk->sub[0].fval;
2272 #else
2273               vals=blk->sub[0].dval;
2274 #endif
2275
2276               if (blk->sub[0].nr != (size_t)nor)
2277                   gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2278               if (bORA || bODA)
2279               {
2280                   for(i=0; i<nor; i++)
2281                       orient[i] += vals[i];
2282               }
2283               if (bODR)
2284               {
2285                   for(i=0; i<nor; i++)
2286                       odrms[i] += sqr(vals[i]-oobs[i]);
2287               }
2288               if (bORT)
2289               {
2290                   fprintf(fort,"  %10f",fr->t);
2291                   for(i=0; i<norsel; i++)
2292                       fprintf(fort," %g",vals[orsel[i]]);
2293                   fprintf(fort,"\n");
2294               }
2295               if (bODT)
2296               {
2297                   fprintf(fodt,"  %10f",fr->t);
2298                   for(i=0; i<norsel; i++)
2299                       fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
2300                   fprintf(fodt,"\n");
2301               }
2302               norfr++;
2303           }
2304          blk = find_block_id_enxframe(fr, enxORT, NULL);
2305          if (bOTEN && blk)
2306          {
2307 #ifndef GMX_DOUBLE
2308              xdr_datatype dt=xdr_datatype_float;
2309 #else
2310              xdr_datatype dt=xdr_datatype_double;
2311 #endif
2312              real *vals;
2313
2314              if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2315                  gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2316 #ifndef GMX_DOUBLE
2317              vals=blk->sub[0].fval;
2318 #else
2319              vals=blk->sub[0].dval;
2320 #endif
2321
2322               if (blk->sub[0].nr != (size_t)(nex*12))
2323                   gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2324                             blk->sub[0].nr/12, nex);
2325               fprintf(foten,"  %10f",fr->t);
2326               for(i=0; i<nex; i++)
2327                   for(j=0; j<(bOvec?12:3); j++)
2328                       fprintf(foten," %g",vals[i*12+j]);
2329               fprintf(foten,"\n");
2330           }
2331         }
2332       }
2333       teller++;
2334     }
2335   } while (bCont && (timecheck == 0));
2336
2337   fprintf(stderr,"\n");
2338   close_enx(fp);
2339   if (out)
2340       ffclose(out);
2341
2342   if (bDRAll)
2343       ffclose(fp_pairs);
2344
2345   if (bORT)
2346       ffclose(fort);
2347   if (bODT)
2348       ffclose(fodt);
2349   if (bORA)
2350   {
2351       out = xvgropen(opt2fn("-ora",NFILE,fnm),
2352                      "Average calculated orientations",
2353                      "Restraint label","",oenv);
2354       if (bOrinst)
2355           fprintf(out,"%s",orinst_sub);
2356       for(i=0; i<nor; i++)
2357           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr);
2358       ffclose(out);
2359   }
2360   if (bODA) {
2361       out = xvgropen(opt2fn("-oda",NFILE,fnm),
2362                      "Average restraint deviation",
2363                      "Restraint label","",oenv);
2364       if (bOrinst)
2365           fprintf(out,"%s",orinst_sub);
2366       for(i=0; i<nor; i++)
2367           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr-oobs[i]);
2368       ffclose(out);
2369   }
2370   if (bODR) {
2371       out = xvgropen(opt2fn("-odr",NFILE,fnm),
2372                      "RMS orientation restraint deviations",
2373                      "Restraint label","",oenv);
2374       if (bOrinst)
2375           fprintf(out,"%s",orinst_sub);
2376       for(i=0; i<nor; i++)
2377           fprintf(out,"%5d  %g\n",or_label[i],sqrt(odrms[i]/norfr));
2378       ffclose(out);
2379   }
2380   if (bOTEN)
2381       ffclose(foten);
2382
2383   if (bDisRe)
2384   {
2385       analyse_disre(opt2fn("-viol",NFILE,fnm),
2386                     teller_disre,violaver,bounds,index,pair,nbounds,oenv);
2387   }
2388   else if (bDHDL)
2389   {
2390       if (fp_dhdl)
2391       {
2392           ffclose(fp_dhdl);
2393           printf("\n\nWrote %d lambda values with %d samples as ",
2394                  dh_lambdas, dh_samples);
2395           if (dh_hists > 0)
2396           {
2397               printf("%d dH histograms ", dh_hists);
2398           }
2399           if (dh_blocks> 0)
2400           {
2401               printf("%d dH data blocks ", dh_blocks);
2402           }
2403           printf("to %s\n", opt2fn("-odh",NFILE,fnm));
2404
2405       }
2406       else
2407       {
2408           gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f",NFILE,fnm));
2409       }
2410
2411   }
2412   else
2413   {
2414       double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2415       analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
2416                    bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
2417                    bVisco,opt2fn("-vis",NFILE,fnm),
2418                    nmol,
2419                    start_step,start_t,frame[cur].step,frame[cur].t,
2420                    time,reftemp,&edat,
2421                    nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
2422                    oenv);
2423       calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
2424                              nbmin,nbmax);
2425   }
2426   if (opt2bSet("-f2",NFILE,fnm)) {
2427       fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
2428           reftemp, nset, set, leg, &edat, time ,oenv);
2429   }
2430
2431   {
2432       const char *nxy = "-nxy";
2433
2434       do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
2435       do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
2436       do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
2437       do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
2438       do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
2439       do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
2440       do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
2441       do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
2442       do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
2443   }
2444   thanx(stderr);
2445
2446   return 0;
2447 }