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, and that the command line option",
1664     "[TT]-fluct_props[tt] is given. The following properties",
1665     "will be computed:[BR]",
1666     "Property                        Energy terms needed[BR]",
1667     "---------------------------------------------------[BR]",
1668     "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp     [BR]",
1669     "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp         [BR]",
1670     "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1671     "Isothermal compressibility:     Vol, Temp          [BR]",
1672     "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1673     "---------------------------------------------------[BR]",
1674     "You always need to set the number of molecules [TT]-nmol[tt].",
1675     "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1676     "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
1677     "When the [TT]-viol[tt] option is set, the time averaged",
1678     "violations are plotted and the running time-averaged and",
1679     "instantaneous sum of violations are recalculated. Additionally",
1680     "running time-averaged and instantaneous distances between",
1681     "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1682
1683     "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1684     "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1685     "The first two options plot the orientation, the last three the",
1686     "deviations of the orientations from the experimental values.",
1687     "The options that end on an 'a' plot the average over time",
1688     "as a function of restraint. The options that end on a 't'",
1689     "prompt the user for restraint label numbers and plot the data",
1690     "as a function of time. Option [TT]-odr[tt] plots the RMS",
1691     "deviation as a function of restraint.",
1692     "When the run used time or ensemble averaged orientation restraints,",
1693     "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1694     "not ensemble-averaged orientations and deviations instead of",
1695     "the time and ensemble averages.[PAR]",
1696
1697     "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1698     "tensor for each orientation restraint experiment. With option",
1699     "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1700
1701     "Option [TT]-odh[tt] extracts and plots the free energy data",
1702     "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1703     "from the [TT]ener.edr[tt] file.[PAR]",
1704
1705     "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1706     "difference with an ideal gas state: [BR]",
1707     "  [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]",
1708     "  [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]",
1709     "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1710     "the average is over the ensemble (or time in a trajectory).",
1711     "Note that this is in principle",
1712     "only correct when averaging over the whole (Boltzmann) ensemble",
1713     "and using the potential energy. This also allows for an entropy",
1714     "estimate using:[BR]",
1715     "  [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]",
1716     "  [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",
1717     "[PAR]",
1718
1719     "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1720     "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1721     "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1722     "files, and the average is over the ensemble A. The running average",
1723     "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1724     "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1725
1726   };
1727   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
1728   static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE,bFluctProps=FALSE;
1729   static int  skip=0,nmol=1,nbmin=5,nbmax=5;
1730   static real reftemp=300.0,ezero=0;
1731   t_pargs pa[] = {
1732     { "-fee",   FALSE, etBOOL,  {&bFee},
1733       "Do a free energy estimate" },
1734     { "-fetemp", FALSE, etREAL,{&reftemp},
1735       "Reference temperature for free energy calculation" },
1736     { "-zero", FALSE, etREAL, {&ezero},
1737       "Subtract a zero-point energy" },
1738     { "-sum",  FALSE, etBOOL, {&bSum},
1739       "Sum the energy terms selected rather than display them all" },
1740     { "-dp",   FALSE, etBOOL, {&bDp},
1741       "Print energies in high precision" },
1742     { "-nbmin", FALSE, etINT, {&nbmin},
1743       "Minimum number of blocks for error estimate" },
1744     { "-nbmax", FALSE, etINT, {&nbmax},
1745       "Maximum number of blocks for error estimate" },
1746     { "-mutot",FALSE, etBOOL, {&bMutot},
1747       "Compute the total dipole moment from the components" },
1748     { "-skip", FALSE, etINT,  {&skip},
1749       "Skip number of frames between data points" },
1750     { "-aver", FALSE, etBOOL, {&bPrAll},
1751       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1752     { "-nmol", FALSE, etINT,  {&nmol},
1753       "Number of molecules in your sample: the energies are divided by this number" },
1754     { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1755       "Compute properties based on energy fluctuations, like heat capacity" },
1756     { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1757       "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1758     { "-fluc", FALSE, etBOOL, {&bFluct},
1759       "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1760     { "-orinst", FALSE, etBOOL, {&bOrinst},
1761       "Analyse instantaneous orientation data" },
1762     { "-ovec", FALSE, etBOOL, {&bOvec},
1763       "Also plot the eigenvectors with [TT]-oten[tt]" }
1764   };
1765   const char* drleg[] = {
1766     "Running average",
1767     "Instantaneous"
1768   };
1769   static const char *setnm[] = {
1770     "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1771     "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1772     "Volume",  "Pressure"
1773   };
1774
1775   FILE       *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1776   FILE       *fp_dhdl=NULL;
1777   FILE       **drout;
1778   ener_file_t fp;
1779   int        timecheck=0;
1780   gmx_mtop_t mtop;
1781   gmx_localtop_t *top=NULL;
1782   t_inputrec ir;
1783   t_energy   **ee;
1784   enerdata_t edat;
1785   gmx_enxnm_t *enm=NULL;
1786   t_enxframe *frame,*fr=NULL;
1787   int        cur=0;
1788 #define NEXT (1-cur)
1789   int        nre,teller,teller_disre,nfr;
1790   gmx_large_int_t start_step;
1791   int        nor=0,nex=0,norfr=0,enx_i=0;
1792   real       start_t;
1793   real       *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1794   int        *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1795   int        nbounds=0,npairs;
1796   gmx_bool       bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN,bDHDL;
1797   gmx_bool       bFoundStart,bCont,bEDR,bVisco;
1798   double     sum,sumaver,sumt,ener,dbl;
1799   double     *time=NULL;
1800   real       Vaver;
1801   int        *set=NULL,i,j,k,nset,sss;
1802   gmx_bool       *bIsEner=NULL;
1803   char       **pairleg,**odtleg,**otenleg;
1804   char       **leg=NULL;
1805   char       **nms;
1806   char       *anm_j,*anm_k,*resnm_j,*resnm_k;
1807   int        resnr_j,resnr_k;
1808   const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1809   char       buf[256];
1810   output_env_t oenv;
1811   t_enxblock *blk=NULL;
1812   t_enxblock *blk_disre=NULL;
1813   int        ndisre=0;
1814   int        dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas=0;
1815
1816   t_filenm   fnm[] = {
1817     { efEDR, "-f",    NULL,      ffREAD  },
1818     { efEDR, "-f2",   NULL,      ffOPTRD },
1819     { efTPX, "-s",    NULL,      ffOPTRD },
1820     { efXVG, "-o",    "energy",  ffWRITE },
1821     { efXVG, "-viol", "violaver",ffOPTWR },
1822     { efXVG, "-pairs","pairs",   ffOPTWR },
1823     { efXVG, "-ora",  "orienta", ffOPTWR },
1824     { efXVG, "-ort",  "orientt", ffOPTWR },
1825     { efXVG, "-oda",  "orideva", ffOPTWR },
1826     { efXVG, "-odr",  "oridevr", ffOPTWR },
1827     { efXVG, "-odt",  "oridevt", ffOPTWR },
1828     { efXVG, "-oten", "oriten",  ffOPTWR },
1829     { efXVG, "-corr", "enecorr", ffOPTWR },
1830     { efXVG, "-vis",  "visco",   ffOPTWR },
1831     { efXVG, "-ravg", "runavgdf",ffOPTWR },
1832     { efXVG, "-odh",  "dhdl"    ,ffOPTWR }
1833   };
1834 #define NFILE asize(fnm)
1835   int     npargs;
1836   t_pargs *ppa;
1837
1838   CopyRight(stderr,argv[0]);
1839   npargs = asize(pa);
1840   ppa    = add_acf_pargs(&npargs,pa);
1841   parse_common_args(&argc,argv,
1842                     PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1843                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1844
1845   bDRAll = opt2bSet("-pairs",NFILE,fnm);
1846   bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1847   bORA   = opt2bSet("-ora",NFILE,fnm);
1848   bORT   = opt2bSet("-ort",NFILE,fnm);
1849   bODA   = opt2bSet("-oda",NFILE,fnm);
1850   bODR   = opt2bSet("-odr",NFILE,fnm);
1851   bODT   = opt2bSet("-odt",NFILE,fnm);
1852   bORIRE = bORA || bORT || bODA || bODR || bODT;
1853   bOTEN  = opt2bSet("-oten",NFILE,fnm);
1854   bDHDL  = opt2bSet("-odh",NFILE,fnm);
1855
1856   nset = 0;
1857
1858   snew(frame,2);
1859   fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1860   do_enxnms(fp,&nre,&enm);
1861
1862   Vaver = -1;
1863
1864   bVisco = opt2bSet("-vis",NFILE,fnm);
1865
1866   if (!bDisRe && !bDHDL)
1867   {
1868       if (bVisco) {
1869           nset=asize(setnm);
1870           snew(set,nset);
1871           /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1872           for(j=0; j<nset; j++) {
1873               for(i=0; i<nre; i++) {
1874                   if (strstr(enm[i].name,setnm[j])) {
1875                       set[j]=i;
1876                       break;
1877                   }
1878               }
1879               if (i == nre) {
1880                   if (gmx_strcasecmp(setnm[j],"Volume")==0) {
1881                       printf("Enter the box volume (" unit_volume "): ");
1882                       if(1 != scanf("%lf",&dbl))
1883                       {
1884                           gmx_fatal(FARGS,"Error reading user input");
1885                       }
1886                       Vaver = dbl;
1887                   } else
1888                       gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1889                                 setnm[j]);
1890               }
1891           }
1892       }
1893       else
1894       {
1895           set=select_by_name(nre,enm,&nset);
1896       }
1897       /* Print all the different units once */
1898       sprintf(buf,"(%s)",enm[set[0]].unit);
1899       for(i=1; i<nset; i++) {
1900           for(j=0; j<i; j++) {
1901               if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1902                   break;
1903               }
1904           }
1905           if (j == i) {
1906               strcat(buf,", (");
1907               strcat(buf,enm[set[i]].unit);
1908               strcat(buf,")");
1909           }
1910       }
1911       out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1912                    oenv);
1913
1914       snew(leg,nset+1);
1915       for(i=0; (i<nset); i++)
1916           leg[i] = enm[set[i]].name;
1917       if (bSum) {
1918           leg[nset]=strdup("Sum");
1919           xvgr_legend(out,nset+1,(const char**)leg,oenv);
1920       }
1921       else
1922           xvgr_legend(out,nset,(const char**)leg,oenv);
1923
1924       snew(bIsEner,nset);
1925       for(i=0; (i<nset); i++) {
1926           bIsEner[i] = FALSE;
1927           for (j=0; (j <= F_ETOT); j++)
1928               bIsEner[i] = bIsEner[i] ||
1929                         (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1930       }
1931
1932       if (bPrAll && nset > 1) {
1933           gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1934       }
1935
1936       time = NULL;
1937
1938       if (bORIRE || bOTEN)
1939           get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1940
1941       if (bORIRE) {
1942           if (bOrinst)
1943               enx_i = enxORI;
1944           else
1945               enx_i = enxOR;
1946
1947           if (bORA || bODA)
1948               snew(orient,nor);
1949           if (bODR)
1950               snew(odrms,nor);
1951           if (bORT || bODT) {
1952               fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1953               fprintf(stderr,"End your selection with 0\n");
1954               j = -1;
1955               orsel = NULL;
1956               do {
1957                   j++;
1958                   srenew(orsel,j+1);
1959                   if(1 != scanf("%d",&(orsel[j])))
1960                   {
1961                       gmx_fatal(FARGS,"Error reading user input");
1962                   }
1963               } while (orsel[j] > 0);
1964               if (orsel[0] == -1) {
1965                   fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1966                   norsel = nor;
1967                   srenew(orsel,nor);
1968                   for(i=0; i<nor; i++)
1969                       orsel[i] = i;
1970               } else {
1971                   /* Build the selection */
1972                   norsel=0;
1973                   for(i=0; i<j; i++) {
1974                       for(k=0; k<nor; k++)
1975                           if (or_label[k] == orsel[i]) {
1976                               orsel[norsel] = k;
1977                               norsel++;
1978                               break;
1979                           }
1980                       if (k == nor)
1981                           fprintf(stderr,"Orientation restraint label %d not found\n",
1982                                   orsel[i]);
1983                   }
1984               }
1985               snew(odtleg,norsel);
1986               for(i=0; i<norsel; i++) {
1987                   snew(odtleg[i],256);
1988                   sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1989               }
1990               if (bORT) {
1991                   fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1992                                 "Time (ps)","",oenv);
1993                   if (bOrinst)
1994                       fprintf(fort,"%s",orinst_sub);
1995                   xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
1996               }
1997               if (bODT) {
1998                   fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1999                                 "Orientation restraint deviation",
2000                                 "Time (ps)","",oenv);
2001                   if (bOrinst)
2002                       fprintf(fodt,"%s",orinst_sub);
2003                   xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
2004               }
2005           }
2006       }
2007       if (bOTEN) {
2008           foten=xvgropen(opt2fn("-oten",NFILE,fnm),
2009                          "Order tensor","Time (ps)","",oenv);
2010           snew(otenleg,bOvec ? nex*12 : nex*3);
2011           for(i=0; i<nex; i++) {
2012               for(j=0; j<3; j++) {
2013                   sprintf(buf,"eig%d",j+1);
2014                   otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
2015               }
2016               if (bOvec) {
2017                   for(j=0; j<9; j++) {
2018                       sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
2019                       otenleg[12*i+3+j] = strdup(buf);
2020                   }
2021               }
2022           }
2023           xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
2024       }
2025   }
2026   else if (bDisRe)
2027   {
2028       nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
2029                          &mtop,&top,&ir);
2030       snew(violaver,npairs);
2031       out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
2032                    "Time (ps)","nm",oenv);
2033       xvgr_legend(out,2,drleg,oenv);
2034       if (bDRAll) {
2035           fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
2036                             "Time (ps)","Distance (nm)",oenv);
2037           if (output_env_get_print_xvgr_codes(oenv))
2038               fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
2039                       ir.dr_tau);
2040       }
2041   }
2042
2043
2044   /* Initiate energies and set them to zero */
2045   edat.nsteps  = 0;
2046   edat.npoints = 0;
2047   edat.nframes = 0;
2048   edat.step    = NULL;
2049   edat.steps   = NULL;
2050   edat.points  = NULL;
2051   snew(edat.s,nset);
2052
2053   /* Initiate counters */
2054   teller       = 0;
2055   teller_disre = 0;
2056   bFoundStart  = FALSE;
2057   start_step   = 0;
2058   start_t      = 0;
2059   do {
2060     /* This loop searches for the first frame (when -b option is given),
2061      * or when this has been found it reads just one energy frame
2062      */
2063     do {
2064       bCont = do_enx(fp,&(frame[NEXT]));
2065
2066       if (bCont) {
2067         timecheck = check_times(frame[NEXT].t);
2068       }
2069     } while (bCont && (timecheck < 0));
2070
2071     if ((timecheck == 0) && bCont) {
2072       /* We read a valid frame, so we can use it */
2073       fr = &(frame[NEXT]);
2074
2075       if (fr->nre > 0) {
2076         /* The frame contains energies, so update cur */
2077         cur  = NEXT;
2078
2079                 if (edat.nframes % 1000 == 0)
2080             {
2081                 srenew(edat.step,edat.nframes+1000);
2082                 memset(&(edat.step[edat.nframes]),0,1000*sizeof(edat.step[0]));
2083                 srenew(edat.steps,edat.nframes+1000);
2084                 memset(&(edat.steps[edat.nframes]),0,1000*sizeof(edat.steps[0]));
2085                 srenew(edat.points,edat.nframes+1000);
2086                 memset(&(edat.points[edat.nframes]),0,1000*sizeof(edat.points[0]));
2087                 for(i=0; i<nset; i++)
2088                 {
2089                     srenew(edat.s[i].ener,edat.nframes+1000);
2090                     memset(&(edat.s[i].ener[edat.nframes]),0,
2091                            1000*sizeof(edat.s[i].ener[0]));
2092
2093                     srenew(edat.s[i].es  ,edat.nframes+1000);
2094                     memset(&(edat.s[i].es[edat.nframes]),0,
2095                            1000*sizeof(edat.s[i].es[0]));
2096                 }
2097             }
2098
2099                 nfr = edat.nframes;
2100             edat.step[nfr] = fr->step;
2101
2102             if (!bFoundStart)
2103             {
2104                 bFoundStart = TRUE;
2105                 /* Initiate the previous step data */
2106                 start_step = fr->step;
2107                 start_t    = fr->t;
2108                 /* Initiate the energy sums */
2109                 edat.steps[nfr]  = 1;
2110                 edat.points[nfr] = 1;
2111                 for(i=0; i<nset; i++)
2112                 {
2113                     sss = set[i];
2114                     edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2115                     edat.s[i].es[nfr].sum2 = 0;
2116                 }
2117                 edat.nsteps  = 1;
2118                 edat.npoints = 1;
2119             }
2120             else
2121             {
2122                 edat.steps[nfr] = fr->nsteps;
2123                 {
2124                     if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2125                     {
2126                         if (fr->nsum <= 1)
2127                         {
2128                             edat.points[nfr] = 1;
2129                             for(i=0; i<nset; i++)
2130                             {
2131                                 sss = set[i];
2132                                 edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2133                                 edat.s[i].es[nfr].sum2 = 0;
2134                             }
2135                             edat.npoints += 1;
2136                         }
2137                         else
2138                         {
2139                             edat.points[nfr] = fr->nsum;
2140                             for(i=0; i<nset; i++)
2141                             {
2142                                 sss = set[i];
2143                                 edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2144                                 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2145                             }
2146                             edat.npoints += fr->nsum;
2147                         }
2148                     }
2149                     else
2150                     {
2151                         /* The interval does not match fr->nsteps:
2152                          * can not do exact averages.
2153                          */
2154                         edat.npoints = 0;
2155                     }
2156                     edat.nsteps = fr->step - start_step + 1;
2157                 }
2158             }
2159             for(i=0; i<nset; i++)
2160             {
2161                 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2162             }
2163       }
2164       /*
2165        * Define distance restraint legends. Can only be done after
2166        * the first frame has been read... (Then we know how many there are)
2167        */
2168       blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
2169       if (bDisRe && bDRAll && !leg && blk_disre)
2170       {
2171           t_iatom   *fa;
2172           t_iparams *ip;
2173
2174           fa = top->idef.il[F_DISRES].iatoms;
2175           ip = top->idef.iparams;
2176           if (blk_disre->nsub != 2 ||
2177               (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2178           {
2179               gmx_incons("Number of disre sub-blocks not equal to 2");
2180           }
2181
2182           ndisre=blk_disre->sub[0].nr ;
2183           if (ndisre != top->idef.il[F_DISRES].nr/3)
2184           {
2185               gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2186                         ndisre,top->idef.il[F_DISRES].nr/3);
2187           }
2188           snew(pairleg,ndisre);
2189           for(i=0; i<ndisre; i++)
2190           {
2191               snew(pairleg[i],30);
2192               j=fa[3*i+1];
2193               k=fa[3*i+2];
2194               gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
2195               gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
2196               sprintf(pairleg[i],"%d %s %d %s (%d)",
2197                       resnr_j,anm_j,resnr_k,anm_k,
2198                       ip[fa[3*i]].disres.label);
2199           }
2200           set=select_it(ndisre,pairleg,&nset);
2201           snew(leg,2*nset);
2202           for(i=0; (i<nset); i++)
2203           {
2204               snew(leg[2*i],32);
2205               sprintf(leg[2*i],  "a %s",pairleg[set[i]]);
2206               snew(leg[2*i+1],32);
2207               sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
2208           }
2209           xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
2210       }
2211
2212       /*
2213        * Store energies for analysis afterwards...
2214        */
2215       if (!bDisRe && !bDHDL && (fr->nre > 0)) {
2216         if (edat.nframes % 1000 == 0) {
2217           srenew(time,edat.nframes+1000);
2218         }
2219         time[edat.nframes] = fr->t;
2220         edat.nframes++;
2221       }
2222       /*
2223        * Printing time, only when we do not want to skip frames
2224        */
2225       if (!skip || teller % skip == 0) {
2226         if (bDisRe) {
2227           /*******************************************
2228            * D I S T A N C E   R E S T R A I N T S
2229            *******************************************/
2230           if (ndisre > 0)
2231           {
2232 #ifndef GMX_DOUBLE
2233             float *disre_rt =     blk_disre->sub[0].fval;
2234             float *disre_rm3tav = blk_disre->sub[1].fval;
2235 #else
2236             double *disre_rt =     blk_disre->sub[0].dval;
2237             double *disre_rm3tav = blk_disre->sub[1].dval;
2238 #endif
2239
2240             print_time(out,fr->t);
2241             if (violaver == NULL)
2242               snew(violaver,ndisre);
2243
2244             /* Subtract bounds from distances, to calculate violations */
2245             calc_violations(disre_rt, disre_rm3tav,
2246                             nbounds,pair,bounds,violaver,&sumt,&sumaver);
2247
2248             fprintf(out,"  %8.4f  %8.4f\n",sumaver,sumt);
2249             if (bDRAll) {
2250               print_time(fp_pairs,fr->t);
2251               for(i=0; (i<nset); i++) {
2252                 sss=set[i];
2253                 fprintf(fp_pairs,"  %8.4f", mypow(disre_rm3tav[sss],minthird));
2254                 fprintf(fp_pairs,"  %8.4f", disre_rt[sss]);
2255               }
2256               fprintf(fp_pairs,"\n");
2257             }
2258             teller_disre++;
2259           }
2260         }
2261         else if (bDHDL)
2262         {
2263             do_dhdl(fr, &fp_dhdl, opt2fn("-odh",NFILE,fnm),
2264                     &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas,
2265                     oenv);
2266         }
2267         /*******************************************
2268          * E N E R G I E S
2269          *******************************************/
2270         else {
2271           if (fr->nre > 0) {
2272             if (bPrAll)
2273             {
2274                 /* We skip frames with single points (usually only the first frame),
2275                  * since they would result in an average plot with outliers.
2276                  */
2277                 if (fr->nsum > 1) {
2278                     print_time(out,fr->t);
2279                      print1(out,bDp,fr->ener[set[0]].e);
2280                      print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
2281                      print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
2282                      fprintf(out,"\n");
2283                 }
2284             }
2285             else
2286             {
2287                 print_time(out,fr->t);
2288                 if (bSum)
2289                 {
2290                     sum = 0;
2291                     for(i=0; i<nset; i++)
2292                     {
2293                         sum += fr->ener[set[i]].e;
2294                     }
2295                     print1(out,bDp,sum/nmol-ezero);
2296                 }
2297                 else
2298                 {
2299                     for(i=0; (i<nset); i++)
2300                     {
2301                         if (bIsEner[i])
2302                         {
2303                             print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
2304                         }
2305                         else
2306                         {
2307                             print1(out,bDp,fr->ener[set[i]].e);
2308                         }
2309                     }
2310                 }
2311                 fprintf(out,"\n");
2312             }
2313           }
2314 #if 0
2315           /* we first count the blocks that have id 0: the orire blocks */
2316           block_orire=0;
2317           for(b=0;b<fr->nblock;b++)
2318           {
2319               if (fr->block[b].id == mde_block_type_orire)
2320                   nblock_orire++;
2321           }
2322 #endif
2323           blk = find_block_id_enxframe(fr, enx_i, NULL);
2324           if (bORIRE && blk)
2325           {
2326 #ifndef GMX_DOUBLE
2327               xdr_datatype dt=xdr_datatype_float;
2328 #else
2329               xdr_datatype dt=xdr_datatype_double;
2330 #endif
2331               real *vals;
2332
2333               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2334               {
2335                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2336               }
2337 #ifndef GMX_DOUBLE
2338               vals=blk->sub[0].fval;
2339 #else
2340               vals=blk->sub[0].dval;
2341 #endif
2342
2343               if (blk->sub[0].nr != (size_t)nor)
2344                   gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2345               if (bORA || bODA)
2346               {
2347                   for(i=0; i<nor; i++)
2348                       orient[i] += vals[i];
2349               }
2350               if (bODR)
2351               {
2352                   for(i=0; i<nor; i++)
2353                       odrms[i] += sqr(vals[i]-oobs[i]);
2354               }
2355               if (bORT)
2356               {
2357                   fprintf(fort,"  %10f",fr->t);
2358                   for(i=0; i<norsel; i++)
2359                       fprintf(fort," %g",vals[orsel[i]]);
2360                   fprintf(fort,"\n");
2361               }
2362               if (bODT)
2363               {
2364                   fprintf(fodt,"  %10f",fr->t);
2365                   for(i=0; i<norsel; i++)
2366                       fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
2367                   fprintf(fodt,"\n");
2368               }
2369               norfr++;
2370           }
2371           blk = find_block_id_enxframe(fr, enxORT, NULL);
2372           if (bOTEN && blk)
2373           {
2374 #ifndef GMX_DOUBLE
2375               xdr_datatype dt=xdr_datatype_float;
2376 #else
2377               xdr_datatype dt=xdr_datatype_double;
2378 #endif
2379               real *vals;
2380
2381               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2382                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2383 #ifndef GMX_DOUBLE
2384               vals=blk->sub[0].fval;
2385 #else
2386               vals=blk->sub[0].dval;
2387 #endif
2388
2389               if (blk->sub[0].nr != (size_t)(nex*12))
2390                   gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2391                             blk->sub[0].nr/12, nex);
2392               fprintf(foten,"  %10f",fr->t);
2393               for(i=0; i<nex; i++)
2394                   for(j=0; j<(bOvec?12:3); j++)
2395                       fprintf(foten," %g",vals[i*12+j]);
2396               fprintf(foten,"\n");
2397           }
2398         }
2399       }
2400       teller++;
2401     }
2402   } while (bCont && (timecheck == 0));
2403
2404   fprintf(stderr,"\n");
2405   close_enx(fp);
2406   if (out)
2407       ffclose(out);
2408
2409   if (bDRAll)
2410       ffclose(fp_pairs);
2411
2412   if (bORT)
2413       ffclose(fort);
2414   if (bODT)
2415       ffclose(fodt);
2416   if (bORA)
2417   {
2418       out = xvgropen(opt2fn("-ora",NFILE,fnm),
2419                      "Average calculated orientations",
2420                      "Restraint label","",oenv);
2421       if (bOrinst)
2422           fprintf(out,"%s",orinst_sub);
2423       for(i=0; i<nor; i++)
2424           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr);
2425       ffclose(out);
2426   }
2427   if (bODA) {
2428       out = xvgropen(opt2fn("-oda",NFILE,fnm),
2429                      "Average restraint deviation",
2430                      "Restraint label","",oenv);
2431       if (bOrinst)
2432           fprintf(out,"%s",orinst_sub);
2433       for(i=0; i<nor; i++)
2434           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr-oobs[i]);
2435       ffclose(out);
2436   }
2437   if (bODR) {
2438       out = xvgropen(opt2fn("-odr",NFILE,fnm),
2439                      "RMS orientation restraint deviations",
2440                      "Restraint label","",oenv);
2441       if (bOrinst)
2442           fprintf(out,"%s",orinst_sub);
2443       for(i=0; i<nor; i++)
2444           fprintf(out,"%5d  %g\n",or_label[i],sqrt(odrms[i]/norfr));
2445       ffclose(out);
2446   }
2447   if (bOTEN)
2448       ffclose(foten);
2449
2450   if (bDisRe)
2451   {
2452       analyse_disre(opt2fn("-viol",NFILE,fnm),
2453                     teller_disre,violaver,bounds,index,pair,nbounds,oenv);
2454   }
2455   else if (bDHDL)
2456   {
2457       if (fp_dhdl)
2458       {
2459           ffclose(fp_dhdl);
2460           printf("\n\nWrote %d lambda values with %d samples as ",
2461                  dh_lambdas, dh_samples);
2462           if (dh_hists > 0)
2463           {
2464               printf("%d dH histograms ", dh_hists);
2465           }
2466           if (dh_blocks> 0)
2467           {
2468               printf("%d dH data blocks ", dh_blocks);
2469           }
2470           printf("to %s\n", opt2fn("-odh",NFILE,fnm));
2471
2472       }
2473       else
2474       {
2475           gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f",NFILE,fnm));
2476       }
2477
2478   }
2479   else
2480   {
2481       double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2482       analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
2483                    bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
2484                    bVisco,opt2fn("-vis",NFILE,fnm),
2485                    nmol,
2486                    start_step,start_t,frame[cur].step,frame[cur].t,
2487                    time,reftemp,&edat,
2488                    nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
2489                    oenv);
2490       if (bFluctProps)
2491           calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
2492                                  nbmin,nbmax);
2493   }
2494   if (opt2bSet("-f2",NFILE,fnm)) {
2495       fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
2496           reftemp, nset, set, leg, &edat, time ,oenv);
2497   }
2498
2499   {
2500       const char *nxy = "-nxy";
2501
2502       do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
2503       do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
2504       do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
2505       do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
2506       do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
2507       do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
2508       do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
2509       do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
2510       do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
2511   }
2512   thanx(stderr);
2513
2514   return 0;
2515 }