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