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