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,
1370                     const char *filename, gmx_bool bDp,
1371                     int *blocks, int *hists, int *samples, int *nlambdas,
1372                     const output_env_t oenv)
1373 {
1374     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
1375     char title[STRLEN],label_x[STRLEN],label_y[STRLEN], legend[STRLEN];
1376     char buf[STRLEN];
1377     gmx_bool first=FALSE;
1378     int nblock_hist=0, nblock_dh=0, nblock_dhcoll=0;
1379     int i,j,k;
1380     /* coll data */
1381     double temp=0, start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
1382     static int setnr=0;
1383     double *native_lambda_vec=NULL;
1384     const char **lambda_components=NULL;
1385     int n_lambda_vec=0;
1386     gmx_bool changing_lambda=FALSE;
1387     int lambda_fep_state;
1388
1389     /* now count the blocks & handle the global dh data */
1390     for(i=0;i<fr->nblock;i++)
1391     {
1392         if (fr->block[i].id == enxDHHIST)
1393         {
1394             nblock_hist++;
1395         }
1396         else if (fr->block[i].id == enxDH)
1397         {
1398             nblock_dh++;
1399         }
1400         else if (fr->block[i].id == enxDHCOLL)
1401         {
1402             nblock_dhcoll++;
1403             if ( (fr->block[i].nsub < 1) ||
1404                  (fr->block[i].sub[0].type != xdr_datatype_double) ||
1405                  (fr->block[i].sub[0].nr < 5))
1406             {
1407                 gmx_fatal(FARGS, "Unexpected block data");
1408             }
1409
1410             /* read the data from the DHCOLL block */
1411             temp =         fr->block[i].sub[0].dval[0];
1412             start_time =   fr->block[i].sub[0].dval[1];
1413             delta_time =   fr->block[i].sub[0].dval[2];
1414             start_lambda = fr->block[i].sub[0].dval[3];
1415             delta_lambda = fr->block[i].sub[0].dval[4];
1416             changing_lambda = (delta_lambda != 0);
1417             if (fr->block[i].nsub > 1)
1418             {
1419                 lambda_fep_state=fr->block[i].sub[1].ival[0];
1420                 if (n_lambda_vec==0)
1421                 {
1422                     n_lambda_vec=fr->block[i].sub[1].ival[1];
1423                 }
1424                 else
1425                 {
1426                     if (n_lambda_vec!=fr->block[i].sub[1].ival[1])
1427                     {
1428                         gmx_fatal(FARGS,
1429                                   "Unexpected change of basis set in lambda");
1430                     }
1431                 }
1432                 if (lambda_components == NULL)
1433                     snew(lambda_components, n_lambda_vec);
1434                 if (native_lambda_vec == NULL)
1435                     snew(native_lambda_vec, n_lambda_vec);
1436                 for(j=0;j<n_lambda_vec;j++)
1437                 {
1438                     native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1439                     lambda_components[j] =
1440                            efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1441                 }
1442             }
1443         }
1444     }
1445
1446     if (nblock_hist == 0 && nblock_dh == 0)
1447     {
1448         /* don't do anything */
1449         return;
1450     }
1451     if (nblock_hist>0 && nblock_dh>0)
1452     {
1453         gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1454     }
1455     if (! *fp_dhdl )
1456     {
1457         if (nblock_dh>0)
1458         {
1459             /* we have standard, non-histogram data --
1460                call open_dhdl to open the file */
1461             /* TODO this is an ugly hack that needs to be fixed: this will only
1462                work if the order of data is always the same and if we're
1463                only using the g_energy compiled with the mdrun that produced
1464                the ener.edr. */
1465             *fp_dhdl=open_dhdl(filename,ir,oenv);
1466         }
1467         else
1468         {
1469             sprintf(title,"N(%s)",deltag);
1470             sprintf(label_x,"%s (%s)",deltag,unit_energy);
1471             sprintf(label_y,"Samples");
1472             *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,oenv);
1473             sprintf(buf,"T = %g (K), %s = %g", temp, lambda, start_lambda);
1474             xvgr_subtitle(*fp_dhdl,buf,oenv);
1475         }
1476     }
1477
1478     (*hists)+=nblock_hist;
1479     (*blocks)+=nblock_dh;
1480     (*nlambdas) = nblock_hist+nblock_dh;
1481
1482     /* write the data */
1483     if (nblock_hist > 0)
1484     {
1485         gmx_large_int_t sum=0;
1486         /* histograms */
1487         for(i=0;i<fr->nblock;i++)
1488         {
1489             t_enxblock *blk=&(fr->block[i]);
1490             if (blk->id==enxDHHIST)
1491             {
1492                 double foreign_lambda, dx;
1493                 gmx_large_int_t x0;
1494                 int nhist, derivative;
1495
1496                 /* check the block types etc. */
1497                 if ( (blk->nsub < 2) ||
1498                      (blk->sub[0].type != xdr_datatype_double) ||
1499                      (blk->sub[1].type != xdr_datatype_large_int) ||
1500                      (blk->sub[0].nr < 2)  ||
1501                      (blk->sub[1].nr < 2) )
1502                 {
1503                     gmx_fatal(FARGS, "Unexpected block data in file");
1504                 }
1505                 foreign_lambda=blk->sub[0].dval[0];
1506                 dx=blk->sub[0].dval[1];
1507                 nhist=blk->sub[1].lval[0];
1508                 derivative=blk->sub[1].lval[1];
1509                 for(j=0;j<nhist;j++)
1510                 {
1511                     const char *lg[1];
1512                     x0=blk->sub[1].lval[2+j];
1513
1514                     if (!derivative)
1515                     {
1516                         sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1517                                 deltag, lambda, foreign_lambda,
1518                                 lambda, start_lambda);
1519                     }
1520                     else
1521                     {
1522                         sprintf(legend, "N(%s | %s=%g)",
1523                                 dhdl, lambda, start_lambda);
1524                     }
1525
1526                     lg[0]=legend;
1527                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1528                     setnr++;
1529                     for(k=0;k<blk->sub[j+2].nr;k++)
1530                     {
1531                         int hist;
1532                         double xmin, xmax;
1533
1534                         hist=blk->sub[j+2].ival[k];
1535                         xmin=(x0+k)*dx;
1536                         xmax=(x0+k+1)*dx;
1537                         fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
1538                                 xmax, hist);
1539                         sum+=hist;
1540                     }
1541                     /* multiple histogram data blocks in one histogram
1542                        mean that the second one is the reverse of the first one:
1543                        for dhdl derivatives, it's important to know both the
1544                        maximum and minimum values */
1545                     dx=-dx;
1546                 }
1547             }
1548         }
1549         (*samples) += (int)(sum/nblock_hist);
1550     }
1551     else
1552     {
1553         /* raw dh */
1554         int len=0;
1555         char **setnames=NULL;
1556         int nnames=nblock_dh;
1557
1558         for(i=0;i<fr->nblock;i++)
1559         {
1560             t_enxblock *blk=&(fr->block[i]);
1561             if (blk->id == enxDH)
1562             {
1563                 if (len == 0)
1564                 {
1565                     len=blk->sub[2].nr;
1566                 }
1567                 else
1568                 {
1569                     if (len!=blk->sub[2].nr)
1570                     {
1571                         gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1572                     }
1573                 }
1574             }
1575         }
1576         (*samples) += len;
1577
1578         for(i=0;i<len;i++)
1579         {
1580             double time=start_time + delta_time*i;
1581
1582             fprintf(*fp_dhdl,"%.4f ", time);
1583
1584             for(j=0;j<fr->nblock;j++)
1585             {
1586                 t_enxblock *blk=&(fr->block[j]);
1587                 if (blk->id == enxDH)
1588                 {
1589                     double value;
1590                     if (blk->sub[2].type == xdr_datatype_float)
1591                     {
1592                         value=blk->sub[2].fval[i];
1593                     }
1594                     else
1595                     {
1596                         value=blk->sub[2].dval[i];
1597                     }
1598                     /* we need to decide which data type it is based on the count*/
1599
1600                     if (j==1 && ir->bExpanded)
1601                     {
1602                         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! */
1603                     } else {
1604                         if (bDp) {
1605                             fprintf(*fp_dhdl," %#.12g", value);   /* print normal precision */
1606                         }
1607                         else
1608                         {
1609                             fprintf(*fp_dhdl," %#.8g", value);   /* print normal precision */
1610                         }
1611                     }
1612                 }
1613             }
1614             fprintf(*fp_dhdl, "\n");
1615         }
1616     }
1617 }
1618
1619
1620 int gmx_energy(int argc,char *argv[])
1621 {
1622   const char *desc[] = {
1623
1624     "[TT]g_energy[tt] extracts energy components or distance restraint",
1625     "data from an energy file. The user is prompted to interactively",
1626     "select the desired energy terms.[PAR]",
1627
1628     "Average, RMSD, and drift are calculated with full precision from the",
1629     "simulation (see printed manual). Drift is calculated by performing",
1630     "a least-squares fit of the data to a straight line. The reported total drift",
1631     "is the difference of the fit at the first and last point.",
1632     "An error estimate of the average is given based on a block averages",
1633     "over 5 blocks using the full-precision averages. The error estimate",
1634     "can be performed over multiple block lengths with the options",
1635     "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1636     "[BB]Note[bb] that in most cases the energy files contains averages over all",
1637     "MD steps, or over many more points than the number of frames in",
1638     "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
1639     "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1640     "file, the statistics mentioned above are simply over the single, per-frame",
1641     "energy values.[PAR]",
1642
1643     "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1644
1645     "Some fluctuation-dependent properties can be calculated provided",
1646     "the correct energy terms are selected, and that the command line option",
1647     "[TT]-fluct_props[tt] is given. The following properties",
1648     "will be computed:[BR]",
1649     "Property                        Energy terms needed[BR]",
1650     "---------------------------------------------------[BR]",
1651     "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp     [BR]",
1652     "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp         [BR]",
1653     "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1654     "Isothermal compressibility:     Vol, Temp          [BR]",
1655     "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1656     "---------------------------------------------------[BR]",
1657     "You always need to set the number of molecules [TT]-nmol[tt].",
1658     "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1659     "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
1660     "When the [TT]-viol[tt] option is set, the time averaged",
1661     "violations are plotted and the running time-averaged and",
1662     "instantaneous sum of violations are recalculated. Additionally",
1663     "running time-averaged and instantaneous distances between",
1664     "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1665
1666     "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1667     "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1668     "The first two options plot the orientation, the last three the",
1669     "deviations of the orientations from the experimental values.",
1670     "The options that end on an 'a' plot the average over time",
1671     "as a function of restraint. The options that end on a 't'",
1672     "prompt the user for restraint label numbers and plot the data",
1673     "as a function of time. Option [TT]-odr[tt] plots the RMS",
1674     "deviation as a function of restraint.",
1675     "When the run used time or ensemble averaged orientation restraints,",
1676     "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1677     "not ensemble-averaged orientations and deviations instead of",
1678     "the time and ensemble averages.[PAR]",
1679
1680     "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1681     "tensor for each orientation restraint experiment. With option",
1682     "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1683
1684     "Option [TT]-odh[tt] extracts and plots the free energy data",
1685     "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1686     "from the [TT]ener.edr[tt] file.[PAR]",
1687
1688     "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1689     "difference with an ideal gas state: [BR]",
1690     "  [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]",
1691     "  [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]",
1692     "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1693     "the average is over the ensemble (or time in a trajectory).",
1694     "Note that this is in principle",
1695     "only correct when averaging over the whole (Boltzmann) ensemble",
1696     "and using the potential energy. This also allows for an entropy",
1697     "estimate using:[BR]",
1698     "  [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]",
1699     "  [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",
1700     "[PAR]",
1701
1702     "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1703     "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1704     "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1705     "files, and the average is over the ensemble A. The running average",
1706     "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1707     "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1708
1709   };
1710   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
1711   static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE,bFluctProps=FALSE;
1712   static int  skip=0,nmol=1,nbmin=5,nbmax=5;
1713   static real reftemp=300.0,ezero=0;
1714   t_pargs pa[] = {
1715     { "-fee",   FALSE, etBOOL,  {&bFee},
1716       "Do a free energy estimate" },
1717     { "-fetemp", FALSE, etREAL,{&reftemp},
1718       "Reference temperature for free energy calculation" },
1719     { "-zero", FALSE, etREAL, {&ezero},
1720       "Subtract a zero-point energy" },
1721     { "-sum",  FALSE, etBOOL, {&bSum},
1722       "Sum the energy terms selected rather than display them all" },
1723     { "-dp",   FALSE, etBOOL, {&bDp},
1724       "Print energies in high precision" },
1725     { "-nbmin", FALSE, etINT, {&nbmin},
1726       "Minimum number of blocks for error estimate" },
1727     { "-nbmax", FALSE, etINT, {&nbmax},
1728       "Maximum number of blocks for error estimate" },
1729     { "-mutot",FALSE, etBOOL, {&bMutot},
1730       "Compute the total dipole moment from the components" },
1731     { "-skip", FALSE, etINT,  {&skip},
1732       "Skip number of frames between data points" },
1733     { "-aver", FALSE, etBOOL, {&bPrAll},
1734       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1735     { "-nmol", FALSE, etINT,  {&nmol},
1736       "Number of molecules in your sample: the energies are divided by this number" },
1737     { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1738       "Compute properties based on energy fluctuations, like heat capacity" },
1739     { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1740       "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1741     { "-fluc", FALSE, etBOOL, {&bFluct},
1742       "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1743     { "-orinst", FALSE, etBOOL, {&bOrinst},
1744       "Analyse instantaneous orientation data" },
1745     { "-ovec", FALSE, etBOOL, {&bOvec},
1746       "Also plot the eigenvectors with [TT]-oten[tt]" }
1747   };
1748   const char* drleg[] = {
1749     "Running average",
1750     "Instantaneous"
1751   };
1752   static const char *setnm[] = {
1753     "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1754     "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1755     "Volume",  "Pressure"
1756   };
1757
1758   FILE       *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1759   FILE       *fp_dhdl=NULL;
1760   FILE       **drout;
1761   ener_file_t fp;
1762   int        timecheck=0;
1763   gmx_mtop_t mtop;
1764   gmx_localtop_t *top=NULL;
1765   t_inputrec ir;
1766   t_energy   **ee;
1767   enerdata_t edat;
1768   gmx_enxnm_t *enm=NULL;
1769   t_enxframe *frame,*fr=NULL;
1770   int        cur=0;
1771 #define NEXT (1-cur)
1772   int        nre,teller,teller_disre,nfr;
1773   gmx_large_int_t start_step;
1774   int        nor=0,nex=0,norfr=0,enx_i=0;
1775   real       start_t;
1776   real       *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1777   int        *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1778   int        nbounds=0,npairs;
1779   gmx_bool       bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN,bDHDL;
1780   gmx_bool       bFoundStart,bCont,bEDR,bVisco;
1781   double     sum,sumaver,sumt,ener,dbl;
1782   double     *time=NULL;
1783   real       Vaver;
1784   int        *set=NULL,i,j,k,nset,sss;
1785   gmx_bool       *bIsEner=NULL;
1786   char       **pairleg,**odtleg,**otenleg;
1787   char       **leg=NULL;
1788   char       **nms;
1789   char       *anm_j,*anm_k,*resnm_j,*resnm_k;
1790   int        resnr_j,resnr_k;
1791   const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1792   char       buf[256];
1793   output_env_t oenv;
1794   t_enxblock *blk=NULL;
1795   t_enxblock *blk_disre=NULL;
1796   int        ndisre=0;
1797   int        dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas=0;
1798
1799   t_filenm   fnm[] = {
1800     { efEDR, "-f",    NULL,      ffREAD  },
1801     { efEDR, "-f2",   NULL,      ffOPTRD },
1802     { efTPX, "-s",    NULL,      ffOPTRD },
1803     { efXVG, "-o",    "energy",  ffWRITE },
1804     { efXVG, "-viol", "violaver",ffOPTWR },
1805     { efXVG, "-pairs","pairs",   ffOPTWR },
1806     { efXVG, "-ora",  "orienta", ffOPTWR },
1807     { efXVG, "-ort",  "orientt", ffOPTWR },
1808     { efXVG, "-oda",  "orideva", ffOPTWR },
1809     { efXVG, "-odr",  "oridevr", ffOPTWR },
1810     { efXVG, "-odt",  "oridevt", ffOPTWR },
1811     { efXVG, "-oten", "oriten",  ffOPTWR },
1812     { efXVG, "-corr", "enecorr", ffOPTWR },
1813     { efXVG, "-vis",  "visco",   ffOPTWR },
1814     { efXVG, "-ravg", "runavgdf",ffOPTWR },
1815     { efXVG, "-odh",  "dhdl"    ,ffOPTWR }
1816   };
1817 #define NFILE asize(fnm)
1818   int     npargs;
1819   t_pargs *ppa;
1820
1821   npargs = asize(pa);
1822   ppa    = add_acf_pargs(&npargs,pa);
1823   parse_common_args(&argc,argv,
1824                     PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1825                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1826
1827   bDRAll = opt2bSet("-pairs",NFILE,fnm);
1828   bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1829   bORA   = opt2bSet("-ora",NFILE,fnm);
1830   bORT   = opt2bSet("-ort",NFILE,fnm);
1831   bODA   = opt2bSet("-oda",NFILE,fnm);
1832   bODR   = opt2bSet("-odr",NFILE,fnm);
1833   bODT   = opt2bSet("-odt",NFILE,fnm);
1834   bORIRE = bORA || bORT || bODA || bODR || bODT;
1835   bOTEN  = opt2bSet("-oten",NFILE,fnm);
1836   bDHDL  = opt2bSet("-odh",NFILE,fnm);
1837
1838   nset = 0;
1839
1840   snew(frame,2);
1841   fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1842   do_enxnms(fp,&nre,&enm);
1843
1844   Vaver = -1;
1845
1846   bVisco = opt2bSet("-vis",NFILE,fnm);
1847
1848   if ((!bDisRe) && (!bDHDL))
1849   {
1850       if (bVisco) {
1851           nset=asize(setnm);
1852           snew(set,nset);
1853           /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1854           for(j=0; j<nset; j++) {
1855               for(i=0; i<nre; i++) {
1856                   if (strstr(enm[i].name,setnm[j])) {
1857                       set[j]=i;
1858                       break;
1859                   }
1860               }
1861               if (i == nre) {
1862                   if (gmx_strcasecmp(setnm[j],"Volume")==0) {
1863                       printf("Enter the box volume (" unit_volume "): ");
1864                       if(1 != scanf("%lf",&dbl))
1865                       {
1866                           gmx_fatal(FARGS,"Error reading user input");
1867                       }
1868                       Vaver = dbl;
1869                   } else
1870                       gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1871                                 setnm[j]);
1872               }
1873           }
1874       }
1875       else
1876       {
1877           set=select_by_name(nre,enm,&nset);
1878       }
1879       /* Print all the different units once */
1880       sprintf(buf,"(%s)",enm[set[0]].unit);
1881       for(i=1; i<nset; i++) {
1882           for(j=0; j<i; j++) {
1883               if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1884                   break;
1885               }
1886           }
1887           if (j == i) {
1888               strcat(buf,", (");
1889               strcat(buf,enm[set[i]].unit);
1890               strcat(buf,")");
1891           }
1892       }
1893       out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1894                    oenv);
1895
1896       snew(leg,nset+1);
1897       for(i=0; (i<nset); i++)
1898           leg[i] = enm[set[i]].name;
1899       if (bSum) {
1900           leg[nset]=strdup("Sum");
1901           xvgr_legend(out,nset+1,(const char**)leg,oenv);
1902       }
1903       else
1904           xvgr_legend(out,nset,(const char**)leg,oenv);
1905
1906       snew(bIsEner,nset);
1907       for(i=0; (i<nset); i++) {
1908           bIsEner[i] = FALSE;
1909           for (j=0; (j <= F_ETOT); j++)
1910               bIsEner[i] = bIsEner[i] ||
1911                         (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1912       }
1913
1914       if (bPrAll && nset > 1) {
1915           gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1916       }
1917
1918       time = NULL;
1919
1920       if (bORIRE || bOTEN)
1921           get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1922
1923       if (bORIRE) {
1924           if (bOrinst)
1925               enx_i = enxORI;
1926           else
1927               enx_i = enxOR;
1928
1929           if (bORA || bODA)
1930               snew(orient,nor);
1931           if (bODR)
1932               snew(odrms,nor);
1933           if (bORT || bODT) {
1934               fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1935               fprintf(stderr,"End your selection with 0\n");
1936               j = -1;
1937               orsel = NULL;
1938               do {
1939                   j++;
1940                   srenew(orsel,j+1);
1941                   if(1 != scanf("%d",&(orsel[j])))
1942                   {
1943                       gmx_fatal(FARGS,"Error reading user input");
1944                   }
1945               } while (orsel[j] > 0);
1946               if (orsel[0] == -1) {
1947                   fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1948                   norsel = nor;
1949                   srenew(orsel,nor);
1950                   for(i=0; i<nor; i++)
1951                       orsel[i] = i;
1952               } else {
1953                   /* Build the selection */
1954                   norsel=0;
1955                   for(i=0; i<j; i++) {
1956                       for(k=0; k<nor; k++)
1957                           if (or_label[k] == orsel[i]) {
1958                               orsel[norsel] = k;
1959                               norsel++;
1960                               break;
1961                           }
1962                       if (k == nor)
1963                           fprintf(stderr,"Orientation restraint label %d not found\n",
1964                                   orsel[i]);
1965                   }
1966               }
1967               snew(odtleg,norsel);
1968               for(i=0; i<norsel; i++) {
1969                   snew(odtleg[i],256);
1970                   sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1971               }
1972               if (bORT) {
1973                   fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1974                                 "Time (ps)","",oenv);
1975                   if (bOrinst)
1976                       fprintf(fort,"%s",orinst_sub);
1977                   xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
1978               }
1979               if (bODT) {
1980                   fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1981                                 "Orientation restraint deviation",
1982                                 "Time (ps)","",oenv);
1983                   if (bOrinst)
1984                       fprintf(fodt,"%s",orinst_sub);
1985                   xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
1986               }
1987           }
1988       }
1989       if (bOTEN) {
1990           foten=xvgropen(opt2fn("-oten",NFILE,fnm),
1991                          "Order tensor","Time (ps)","",oenv);
1992           snew(otenleg,bOvec ? nex*12 : nex*3);
1993           for(i=0; i<nex; i++) {
1994               for(j=0; j<3; j++) {
1995                   sprintf(buf,"eig%d",j+1);
1996                   otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
1997               }
1998               if (bOvec) {
1999                   for(j=0; j<9; j++) {
2000                       sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
2001                       otenleg[12*i+3+j] = strdup(buf);
2002                   }
2003               }
2004           }
2005           xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
2006       }
2007   }
2008   else if (bDisRe)
2009   {
2010       nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
2011                          &mtop,&top,&ir);
2012       snew(violaver,npairs);
2013       out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
2014                    "Time (ps)","nm",oenv);
2015       xvgr_legend(out,2,drleg,oenv);
2016       if (bDRAll) {
2017           fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
2018                             "Time (ps)","Distance (nm)",oenv);
2019           if (output_env_get_print_xvgr_codes(oenv))
2020               fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2021                       ir.dr_tau);
2022       }
2023   } else if (bDHDL) {
2024       get_dhdl_parms(ftp2fn(efTPX,NFILE,fnm),&ir);
2025   }
2026
2027    /* Initiate energies and set them to zero */
2028    edat.nsteps  = 0;
2029    edat.npoints = 0;
2030    edat.nframes = 0;
2031    edat.step    = NULL;
2032    edat.steps   = NULL;
2033    edat.points  = NULL;
2034    snew(edat.s,nset);
2035
2036    /* Initiate counters */
2037    teller       = 0;
2038    teller_disre = 0;
2039    bFoundStart  = FALSE;
2040    start_step   = 0;
2041    start_t      = 0;
2042    do {
2043      /* This loop searches for the first frame (when -b option is given),
2044       * or when this has been found it reads just one energy frame
2045       */
2046      do {
2047          bCont = do_enx(fp,&(frame[NEXT]));
2048          if (bCont) {
2049              timecheck = check_times(frame[NEXT].t);
2050          }
2051      } while (bCont && (timecheck < 0));
2052
2053      if ((timecheck == 0) && bCont) {
2054        /* We read a valid frame, so we can use it */
2055        fr = &(frame[NEXT]);
2056
2057        if (fr->nre > 0) {
2058          /* The frame contains energies, so update cur */
2059          cur  = NEXT;
2060
2061              if (edat.nframes % 1000 == 0)
2062              {
2063                  srenew(edat.step,edat.nframes+1000);
2064                  memset(&(edat.step[edat.nframes]),0,1000*sizeof(edat.step[0]));
2065                  srenew(edat.steps,edat.nframes+1000);
2066                  memset(&(edat.steps[edat.nframes]),0,1000*sizeof(edat.steps[0]));
2067                  srenew(edat.points,edat.nframes+1000);
2068                  memset(&(edat.points[edat.nframes]),0,1000*sizeof(edat.points[0]));
2069
2070                  for(i=0; i<nset; i++)
2071                  {
2072                      srenew(edat.s[i].ener,edat.nframes+1000);
2073                      memset(&(edat.s[i].ener[edat.nframes]),0,
2074                             1000*sizeof(edat.s[i].ener[0]));
2075                      srenew(edat.s[i].es  ,edat.nframes+1000);
2076                      memset(&(edat.s[i].es[edat.nframes]),0,
2077                             1000*sizeof(edat.s[i].es[0]));
2078                  }
2079              }
2080
2081              nfr = edat.nframes;
2082              edat.step[nfr] = fr->step;
2083
2084              if (!bFoundStart)
2085              {
2086                  bFoundStart = TRUE;
2087                  /* Initiate the previous step data */
2088                  start_step = fr->step;
2089                  start_t    = fr->t;
2090                  /* Initiate the energy sums */
2091                  edat.steps[nfr]  = 1;
2092                  edat.points[nfr] = 1;
2093                  for(i=0; i<nset; i++)
2094                  {
2095                      sss = set[i];
2096                      edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2097                      edat.s[i].es[nfr].sum2 = 0;
2098                  }
2099                  edat.nsteps  = 1;
2100                  edat.npoints = 1;
2101              }
2102              else
2103              {
2104                  edat.steps[nfr] = fr->nsteps;
2105                  {
2106                      if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2107                      {
2108                          if (fr->nsum <= 1)
2109                          {
2110                              edat.points[nfr] = 1;
2111                              for(i=0; i<nset; i++)
2112                              {
2113                                  sss = set[i];
2114                                  edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2115                                  edat.s[i].es[nfr].sum2 = 0;
2116                              }
2117                              edat.npoints += 1;
2118                          }
2119                          else
2120                          {
2121                              edat.points[nfr] = fr->nsum;
2122                              for(i=0; i<nset; i++)
2123                              {
2124                                  sss = set[i];
2125                                  edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2126                                  edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2127                              }
2128                              edat.npoints += fr->nsum;
2129                          }
2130                      }
2131                      else
2132                      {
2133                          /* The interval does not match fr->nsteps:
2134                           * can not do exact averages.
2135                           */
2136                          edat.npoints = 0;
2137                      }
2138                      edat.nsteps = fr->step - start_step + 1;
2139                  }
2140              }
2141              for(i=0; i<nset; i++)
2142              {
2143                  edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2144              }
2145        }
2146        /*
2147         * Define distance restraint legends. Can only be done after
2148         * the first frame has been read... (Then we know how many there are)
2149         */
2150        blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
2151        if (bDisRe && bDRAll && !leg && blk_disre)
2152        {
2153            t_iatom   *fa;
2154            t_iparams *ip;
2155
2156            fa = top->idef.il[F_DISRES].iatoms;
2157            ip = top->idef.iparams;
2158            if (blk_disre->nsub != 2 ||
2159                (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2160            {
2161                gmx_incons("Number of disre sub-blocks not equal to 2");
2162            }
2163
2164            ndisre=blk_disre->sub[0].nr ;
2165            if (ndisre != top->idef.il[F_DISRES].nr/3)
2166            {
2167                gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2168                          ndisre,top->idef.il[F_DISRES].nr/3);
2169            }
2170            snew(pairleg,ndisre);
2171            for(i=0; i<ndisre; i++)
2172            {
2173                snew(pairleg[i],30);
2174                j=fa[3*i+1];
2175                k=fa[3*i+2];
2176                gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
2177                gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
2178                sprintf(pairleg[i],"%d %s %d %s (%d)",
2179                        resnr_j,anm_j,resnr_k,anm_k,
2180                        ip[fa[3*i]].disres.label);
2181            }
2182            set=select_it(ndisre,pairleg,&nset);
2183            snew(leg,2*nset);
2184            for(i=0; (i<nset); i++)
2185            {
2186                snew(leg[2*i],32);
2187                sprintf(leg[2*i],  "a %s",pairleg[set[i]]);
2188                snew(leg[2*i+1],32);
2189                sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
2190            }
2191            xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
2192        }
2193
2194        /*
2195         * Store energies for analysis afterwards...
2196         */
2197        if (!bDisRe && !bDHDL && (fr->nre > 0)) {
2198            if (edat.nframes % 1000 == 0) {
2199                srenew(time,edat.nframes+1000);
2200            }
2201            time[edat.nframes] = fr->t;
2202            edat.nframes++;
2203        }
2204        /*
2205         * Printing time, only when we do not want to skip frames
2206         */
2207        if (!skip || teller % skip == 0) {
2208      if (bDisRe) {
2209        /*******************************************
2210         * D I S T A N C E   R E S T R A I N T S
2211         *******************************************/
2212        if (ndisre > 0)
2213            {
2214  #ifndef GMX_DOUBLE
2215              float *disre_rt =     blk_disre->sub[0].fval;
2216              float *disre_rm3tav = blk_disre->sub[1].fval;
2217  #else
2218              double *disre_rt =     blk_disre->sub[0].dval;
2219              double *disre_rm3tav = blk_disre->sub[1].dval;
2220  #endif
2221
2222          print_time(out,fr->t);
2223          if (violaver == NULL)
2224            snew(violaver,ndisre);
2225
2226          /* Subtract bounds from distances, to calculate violations */
2227          calc_violations(disre_rt, disre_rm3tav,
2228                  nbounds,pair,bounds,violaver,&sumt,&sumaver);
2229
2230          fprintf(out,"  %8.4f  %8.4f\n",sumaver,sumt);
2231          if (bDRAll) {
2232            print_time(fp_pairs,fr->t);
2233            for(i=0; (i<nset); i++) {
2234          sss=set[i];
2235          fprintf(fp_pairs,"  %8.4f", mypow(disre_rm3tav[sss],minthird));
2236          fprintf(fp_pairs,"  %8.4f", disre_rt[sss]);
2237            }
2238            fprintf(fp_pairs,"\n");
2239          }
2240          teller_disre++;
2241        }
2242      }
2243      else if (bDHDL)
2244      {
2245          do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh",NFILE,fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2246      }
2247
2248      /*******************************************
2249       * E N E R G I E S
2250       *******************************************/
2251      else {
2252          if (fr->nre > 0) {
2253              if (bPrAll)
2254              {
2255                  /* We skip frames with single points (usually only the first frame),
2256                   * since they would result in an average plot with outliers.
2257                   */
2258                  if (fr->nsum > 1) {
2259                      print_time(out,fr->t);
2260                       print1(out,bDp,fr->ener[set[0]].e);
2261                       print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
2262                       print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
2263                       fprintf(out,"\n");
2264                  }
2265              }
2266              else
2267              {
2268                  print_time(out,fr->t);
2269                  if (bSum)
2270                  {
2271                      sum = 0;
2272                      for(i=0; i<nset; i++)
2273                      {
2274                          sum += fr->ener[set[i]].e;
2275                      }
2276                      print1(out,bDp,sum/nmol-ezero);
2277                  }
2278                  else
2279                  {
2280                      for(i=0; (i<nset); i++)
2281                      {
2282                          if (bIsEner[i])
2283                          {
2284                              print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
2285                          }
2286                          else
2287                          {
2288                              print1(out,bDp,fr->ener[set[i]].e);
2289                          }
2290                      }
2291                  }
2292                  fprintf(out,"\n");
2293              }
2294          }
2295           blk = find_block_id_enxframe(fr, enx_i, NULL);
2296           if (bORIRE && blk)
2297           {
2298 #ifndef GMX_DOUBLE
2299               xdr_datatype dt=xdr_datatype_float;
2300 #else
2301               xdr_datatype dt=xdr_datatype_double;
2302 #endif
2303               real *vals;
2304
2305               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2306               {
2307                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2308               }
2309 #ifndef GMX_DOUBLE
2310               vals=blk->sub[0].fval;
2311 #else
2312               vals=blk->sub[0].dval;
2313 #endif
2314
2315               if (blk->sub[0].nr != (size_t)nor)
2316                   gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2317               if (bORA || bODA)
2318               {
2319                   for(i=0; i<nor; i++)
2320                       orient[i] += vals[i];
2321               }
2322               if (bODR)
2323               {
2324                   for(i=0; i<nor; i++)
2325                       odrms[i] += sqr(vals[i]-oobs[i]);
2326               }
2327               if (bORT)
2328               {
2329                   fprintf(fort,"  %10f",fr->t);
2330                   for(i=0; i<norsel; i++)
2331                       fprintf(fort," %g",vals[orsel[i]]);
2332                   fprintf(fort,"\n");
2333               }
2334               if (bODT)
2335               {
2336                   fprintf(fodt,"  %10f",fr->t);
2337                   for(i=0; i<norsel; i++)
2338                       fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
2339                   fprintf(fodt,"\n");
2340               }
2341               norfr++;
2342           }
2343          blk = find_block_id_enxframe(fr, enxORT, NULL);
2344          if (bOTEN && blk)
2345          {
2346 #ifndef GMX_DOUBLE
2347              xdr_datatype dt=xdr_datatype_float;
2348 #else
2349              xdr_datatype dt=xdr_datatype_double;
2350 #endif
2351              real *vals;
2352
2353              if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2354                  gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2355 #ifndef GMX_DOUBLE
2356              vals=blk->sub[0].fval;
2357 #else
2358              vals=blk->sub[0].dval;
2359 #endif
2360
2361               if (blk->sub[0].nr != (size_t)(nex*12))
2362                   gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2363                             blk->sub[0].nr/12, nex);
2364               fprintf(foten,"  %10f",fr->t);
2365               for(i=0; i<nex; i++)
2366                   for(j=0; j<(bOvec?12:3); j++)
2367                       fprintf(foten," %g",vals[i*12+j]);
2368               fprintf(foten,"\n");
2369           }
2370         }
2371       }
2372       teller++;
2373     }
2374   } while (bCont && (timecheck == 0));
2375
2376   fprintf(stderr,"\n");
2377   close_enx(fp);
2378   if (out)
2379       ffclose(out);
2380
2381   if (bDRAll)
2382       ffclose(fp_pairs);
2383
2384   if (bORT)
2385       ffclose(fort);
2386   if (bODT)
2387       ffclose(fodt);
2388   if (bORA)
2389   {
2390       out = xvgropen(opt2fn("-ora",NFILE,fnm),
2391                      "Average calculated orientations",
2392                      "Restraint label","",oenv);
2393       if (bOrinst)
2394           fprintf(out,"%s",orinst_sub);
2395       for(i=0; i<nor; i++)
2396           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr);
2397       ffclose(out);
2398   }
2399   if (bODA) {
2400       out = xvgropen(opt2fn("-oda",NFILE,fnm),
2401                      "Average restraint deviation",
2402                      "Restraint label","",oenv);
2403       if (bOrinst)
2404           fprintf(out,"%s",orinst_sub);
2405       for(i=0; i<nor; i++)
2406           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr-oobs[i]);
2407       ffclose(out);
2408   }
2409   if (bODR) {
2410       out = xvgropen(opt2fn("-odr",NFILE,fnm),
2411                      "RMS orientation restraint deviations",
2412                      "Restraint label","",oenv);
2413       if (bOrinst)
2414           fprintf(out,"%s",orinst_sub);
2415       for(i=0; i<nor; i++)
2416           fprintf(out,"%5d  %g\n",or_label[i],sqrt(odrms[i]/norfr));
2417       ffclose(out);
2418   }
2419   if (bOTEN)
2420       ffclose(foten);
2421
2422   if (bDisRe)
2423   {
2424       analyse_disre(opt2fn("-viol",NFILE,fnm),
2425                     teller_disre,violaver,bounds,index,pair,nbounds,oenv);
2426   }
2427   else if (bDHDL)
2428   {
2429       if (fp_dhdl)
2430       {
2431           ffclose(fp_dhdl);
2432           printf("\n\nWrote %d lambda values with %d samples as ",
2433                  dh_lambdas, dh_samples);
2434           if (dh_hists > 0)
2435           {
2436               printf("%d dH histograms ", dh_hists);
2437           }
2438           if (dh_blocks> 0)
2439           {
2440               printf("%d dH data blocks ", dh_blocks);
2441           }
2442           printf("to %s\n", opt2fn("-odh",NFILE,fnm));
2443
2444       }
2445       else
2446       {
2447           gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f",NFILE,fnm));
2448       }
2449
2450   }
2451   else
2452   {
2453       double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2454       analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
2455                    bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
2456                    bVisco,opt2fn("-vis",NFILE,fnm),
2457                    nmol,
2458                    start_step,start_t,frame[cur].step,frame[cur].t,
2459                    time,reftemp,&edat,
2460                    nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
2461                    oenv);
2462       if (bFluctProps)
2463           calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
2464                                  nbmin,nbmax);
2465   }
2466   if (opt2bSet("-f2",NFILE,fnm)) {
2467       fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
2468           reftemp, nset, set, leg, &edat, time ,oenv);
2469   }
2470
2471   {
2472       const char *nxy = "-nxy";
2473
2474       do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
2475       do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
2476       do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
2477       do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
2478       do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
2479       do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
2480       do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
2481       do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
2482       do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
2483   }
2484   thanx(stderr);
2485
2486   return 0;
2487 }