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