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 analyse_ener(gmx_bool bCorr,const char *corrfn,
 844                          gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,gmx_bool bTempFluct,
 845                          gmx_bool bVisco,const char *visfn,int nmol,
 846                          int nconstr,
 847                          gmx_large_int_t start_step,double start_t,
 848                          gmx_large_int_t step,double t,
 849                          double time[], real reftemp,
 850                          enerdata_t *edat,
 851                          int nset,int set[],gmx_bool *bIsEner,
 852                          char **leg,gmx_enxnm_t *enm,
 853                          real Vaver,real ezero,
 854                          int nbmin,int nbmax,
 855                          const output_env_t oenv)
 856 {
 857   FILE *fp;
 858   /* Check out the printed manual for equations! */
 859   double Dt,aver,stddev,errest,delta_t,totaldrift;
 860   enerdata_t *esum=NULL;
 861   real xxx,integral,intBulk;
 862   real sfrac,oldfrac,diffsum,diffav,fstep,pr_aver,pr_stddev,pr_errest;
 863   double beta=0,expE,expEtot,*fee=NULL;
 864   gmx_large_int_t nsteps;
 865   int  nexact,nnotexact;
 866   double x1m,x1mk;
 867   real Temp=-1,Pres=-1,VarV=-1,VarT=-1,VarEtot=-1,AvEtot=0,VarEnthalpy=-1;
 868   int  i,j,nout;
 869   real chi2;
 870   char buf[256],eebuf[100];
 871
 872   nsteps  = step - start_step + 1;
 873   if (nsteps < 1) {
 874     fprintf(stdout,"Not enough steps (%s) for statistics\n",
 875             gmx_step_str(nsteps,buf));
 876   }
 877   else {
 878     /* Calculate the time difference */
 879     delta_t = t - start_t;
 880
 881     fprintf(stdout,"\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
 882             gmx_step_str(nsteps,buf),start_t,t,nset);
 883
 884     calc_averages(nset,edat,nbmin,nbmax);
 885
 886     if (bSum) {
 887         esum = calc_sum(nset,edat,nbmin,nbmax);
 888     }
 889
 890     if (edat->npoints == 0) {
 891       nexact    = 0;
 892       nnotexact = nset;
 893     } else {
 894       nexact    = 0;
 895       nnotexact = 0;
 896       for(i=0; (i<nset); i++) {
 897         if (edat->s[i].bExactStat) {
 898           nexact++;
 899         } else {
 900           nnotexact++;
 901         }
 902       }
 903     }
 904
 905     if (nnotexact == 0) {
 906       fprintf(stdout,"All statistics are over %s points\n",
 907               gmx_step_str(edat->npoints,buf));
 908     } else if (nexact == 0 || edat->npoints == edat->nframes) {
 909       fprintf(stdout,"All statistics are over %d points (frames)\n",
 910               edat->nframes);
 911     } else {
 912       fprintf(stdout,"The term%s",nnotexact==1 ? "" : "s");
 913       for(i=0; (i<nset); i++) {
 914         if (!edat->s[i].bExactStat) {
 915           fprintf(stdout," '%s'",leg[i]);
 916         }
 917       }
 918       fprintf(stdout," %s has statistics over %d points (frames)\n",
 919               nnotexact==1 ? "is" : "are",edat->nframes);
 920       fprintf(stdout,"All other statistics are over %s points\n",
 921               gmx_step_str(edat->npoints,buf));
 922     }
 923     fprintf(stdout,"\n");
 924
 925     fprintf(stdout,"%-24s %10s %10s %10s %10s",
 926             "Energy","Average","Err.Est.","RMSD","Tot-Drift");
 927     if (bFee)
 928       fprintf(stdout,"  %10s\n","-kT ln<e^(E/kT)>");
 929     else
 930       fprintf(stdout,"\n");
 931     fprintf(stdout,"-------------------------------------------------------------------------------\n");
 932
 933     /* Initiate locals, only used with -sum */
 934     expEtot=0;
 935     if (bFee) {
 936       beta = 1.0/(BOLTZ*reftemp);
 937       snew(fee,nset);
 938     }
 939     for(i=0; (i<nset); i++) {
 940       aver   = edat->s[i].av;
 941       stddev = edat->s[i].rmsd;
 942       errest = edat->s[i].ee;
 943
 944       if (bFee) {
 945         expE = 0;
 946         for(j=0; (j<edat->nframes); j++) {
 947           expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
 948         }
 949         if (bSum)
 950           expEtot+=expE/edat->nframes;
 951
 952         fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
 953       }
 954       if (strstr(leg[i],"empera") != NULL) {
 955         VarT = sqr(stddev);
 956         Temp = aver;
 957       } else if (strstr(leg[i],"olum") != NULL) {
 958         VarV = sqr(stddev);
 959         Vaver= aver;
 960       } else if (strstr(leg[i],"essure") != NULL) {
 961         Pres = aver;
 962       } else if (strstr(leg[i],"otal") != NULL) {
 963         VarEtot = sqr(stddev);
 964         AvEtot = aver;
 965       } else if (strstr(leg[i],"nthalpy") != NULL) {
 966         VarEnthalpy = sqr(stddev);
 967       }
 968       if (bIsEner[i]) {
 969         pr_aver   = aver/nmol-ezero;
 970         pr_stddev = stddev/nmol;
 971         pr_errest = errest/nmol;
 972       }
 973       else {
 974         pr_aver   = aver;
 975         pr_stddev = stddev;
 976         pr_errest = errest;
 977       }
 978
 979       /* Multiply the slope in steps with the number of steps taken */
 980       totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
 981       if (bIsEner[i])
 982       {
 983           totaldrift /= nmol;
 984       }
 985
 986       fprintf(stdout,"%-24s %10g %10s %10g %10g",
 987               leg[i],pr_aver,ee_pr(pr_errest,eebuf),pr_stddev,totaldrift);
 988       if (bFee)
 989         fprintf(stdout,"  %10g",fee[i]);
 990
 991       fprintf(stdout,"  (%s)\n",enm[set[i]].unit);
 992
 993       if (bFluct) {
 994         for(j=0; (j<edat->nframes); j++)
 995           edat->s[i].ener[j] -= aver;
 996       }
 997     }
 998     if (bSum) {
 999         totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1000       fprintf(stdout,"%-24s %10g %10s %10s %10g  (%s)",
1001               "Total",esum->s[0].av/nmol,ee_pr(esum->s[0].ee/nmol,eebuf),
1002               "--",totaldrift/nmol,enm[set[0]].unit);
1003       /* pr_aver,pr_stddev,a,totaldrift */
1004       if (bFee)
1005         fprintf(stdout,"  %10g  %10g\n",
1006                 log(expEtot)/beta + esum->s[0].av/nmol,log(expEtot)/beta);
1007       else
1008         fprintf(stdout,"\n");
1009     }
1010     if (bTempFluct && Temp != -1) {
1011       printf("\nTemperature dependent fluctuation properties at T = %g. #constr/mol = %d\n",Temp,nconstr);
1012       if (nmol < 2)
1013         printf("Warning: nmol = %d, this may not be what you want.\n",
1014                nmol);
1015       if (VarV != -1) {
1016         real tmp = VarV/(Vaver*BOLTZ*Temp*PRESFAC);
1017
1018         printf("Isothermal Compressibility: %10g /%s\n",
1019                tmp,unit_pres_bar);
1020         printf("Adiabatic bulk modulus:     %10g  %s\n",
1021                1.0/tmp,unit_pres_bar);
1022       }
1023       if (VarEnthalpy != -1) {
1024         real Cp = 1000*((VarEnthalpy/nmol)/(BOLTZ*Temp*Temp) -
1025                         0.5*BOLTZ*nconstr);
1026         printf("Heat capacity at constant pressure Cp: %10g J/mol K\n",Cp);
1027       }
1028       if ((VarV != -1) && (VarEnthalpy != -1)) {
1029         real aP = (sqrt(VarEnthalpy*VarV/nmol))/(BOLTZ*Vaver*Temp*Temp);
1030         printf("Thermal expansion coefficient alphaP: %10g 1/K\n",aP);
1031       }
1032       if ((VarV == -1) && (VarEtot != -1)) {
1033         real Cv = 1000*((VarEtot/nmol)/(BOLTZ*Temp*Temp) -
1034                         0.5*BOLTZ*nconstr);
1035         printf("Heat capacity at constant volume Cv: %10g J/mol K\n",Cv);
1036       }
1037       please_cite(stdout,"Allen1987a");
1038     }
1039     /* Do correlation function */
1040     if (edat->nframes > 1)
1041     {
1042         Dt = delta_t/(edat->nframes - 1);
1043     }
1044     else
1045     {
1046         Dt = 0;
1047     }
1048     if (bVisco) {
1049       const char* leg[] = { "Shear", "Bulk" };
1050       real factor;
1051       real **eneset;
1052       real **enesum;
1053
1054       /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1055
1056       /* Symmetrise tensor! (and store in first three elements)
1057        * And subtract average pressure!
1058        */
1059       snew(eneset,12);
1060       for(i=0; i<12; i++) {
1061           snew(eneset[i],edat->nframes);
1062       }
1063       snew(enesum,3);
1064       for(i=0; i<3; i++) {
1065         snew(enesum[i],edat->nframes);
1066       }
1067       for(i=0; (i<edat->nframes); i++) {
1068         eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1069         eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1070         eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1071         for(j=3; j<=11; j++) {
1072           eneset[j][i] = edat->s[j].ener[i];
1073         }
1074         eneset[11][i] -= Pres;
1075         enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1076         enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1077         enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1078       }
1079
1080       einstein_visco("evisco.xvg","eviscoi.xvg",
1081                      3,edat->nframes,enesum,Vaver,Temp,nsteps,time,oenv);
1082
1083       /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1084       /* Do it for shear viscosity */
1085       strcpy(buf,"Shear Viscosity");
1086       low_do_autocorr(corrfn,oenv,buf,edat->nframes,3,
1087                       (edat->nframes+1)/2,eneset,Dt,
1088                       eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1089
1090       /* Now for bulk viscosity */
1091       strcpy(buf,"Bulk Viscosity");
1092       low_do_autocorr(corrfn,oenv,buf,edat->nframes,1,
1093                       (edat->nframes+1)/2,&(eneset[11]),Dt,
1094                       eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1095
1096       factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1097       fp=xvgropen(visfn,buf,"Time (ps)","\\8h\\4 (cp)",oenv);
1098       xvgr_legend(fp,asize(leg),leg,oenv);
1099
1100       /* Use trapezium rule for integration */
1101       integral = 0;
1102       intBulk  = 0;
1103       nout = get_acfnout();
1104       if ((nout < 2) || (nout >= edat->nframes/2))
1105           nout = edat->nframes/2;
1106       for(i=1; (i<nout); i++)
1107       {
1108           integral += 0.5*(eneset[0][i-1]  + eneset[0][i])*factor;
1109           intBulk  += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1110           fprintf(fp,"%10g  %10g  %10g\n",(i*Dt),integral,intBulk);
1111       }
1112       ffclose(fp);
1113     }
1114     else if (bCorr) {
1115       if (bFluct)
1116           strcpy(buf,"Autocorrelation of Energy Fluctuations");
1117       else
1118           strcpy(buf,"Energy Autocorrelation");
1119 #if 0
1120       do_autocorr(corrfn,oenv,buf,edat->nframes,
1121                   bSum ? 1                 : nset,
1122                   bSum ? &edat->s[nset-1].ener : eneset,
1123                   (delta_t/edat->nframes),eacNormal,FALSE);
1124 #endif
1125     }
1126   }
1127 }
1128
1129 static void print_time(FILE *fp,double t)
1130 {
1131   fprintf(fp,"%12.6f",t);
1132 }
1133
1134 static void print1(FILE *fp,gmx_bool bDp,real e)
1135 {
1136   if (bDp)
1137     fprintf(fp,"  %16.12f",e);
1138   else
1139     fprintf(fp,"  %10.6f",e);
1140 }
1141
1142 static void fec(const char *ene2fn, const char *runavgfn,
1143                 real reftemp, int nset, int set[], char *leg[],
1144                 enerdata_t *edat, double time[],
1145                 const output_env_t oenv)
1146 {
1147   const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1148                            "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1149   FILE *fp;
1150   ener_file_t enx;
1151   int  nre,timecheck,step,nenergy,nenergy2,maxenergy;
1152   int  i,j;
1153   gmx_bool bCont;
1154   real aver, beta;
1155   real **eneset2;
1156   double dE, sum;
1157   gmx_enxnm_t *enm=NULL;
1158   t_enxframe *fr;
1159   char buf[22];
1160
1161   /* read second energy file */
1162   snew(fr,1);
1163   enm = NULL;
1164   enx = open_enx(ene2fn,"r");
1165   do_enxnms(enx,&(fr->nre),&enm);
1166
1167   snew(eneset2,nset+1);
1168   nenergy2=0;
1169   maxenergy=0;
1170   timecheck=0;
1171   do {
1172     /* This loop searches for the first frame (when -b option is given),
1173      * or when this has been found it reads just one energy frame
1174      */
1175     do {
1176       bCont = do_enx(enx,fr);
1177
1178       if (bCont)
1179         timecheck = check_times(fr->t);
1180
1181     } while (bCont && (timecheck < 0));
1182
1183     /* Store energies for analysis afterwards... */
1184     if ((timecheck == 0) && bCont) {
1185       if (fr->nre > 0) {
1186         if ( nenergy2 >= maxenergy ) {
1187           maxenergy += 1000;
1188           for(i=0; i<=nset; i++)
1189             srenew(eneset2[i],maxenergy);
1190         }
1191         if (fr->t != time[nenergy2])
1192           fprintf(stderr,"\nWARNING time mismatch %g!=%g at frame %s\n",
1193                   fr->t, time[nenergy2], gmx_step_str(fr->step,buf));
1194         for(i=0; i<nset; i++)
1195           eneset2[i][nenergy2] = fr->ener[set[i]].e;
1196         nenergy2++;
1197       }
1198     }
1199   } while (bCont && (timecheck == 0));
1200
1201   /* check */
1202   if (edat->nframes != nenergy2) {
1203     fprintf(stderr,"\nWARNING file length mismatch %d!=%d\n",
1204             edat->nframes,nenergy2);
1205   }
1206   nenergy = min(edat->nframes,nenergy2);
1207
1208   /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1209   fp=NULL;
1210   if (runavgfn) {
1211     fp=xvgropen(runavgfn,"Running average free energy difference",
1212                 "Time (" unit_time ")","\\8D\\4E (" unit_energy ")",oenv);
1213     xvgr_legend(fp,asize(ravgleg),ravgleg,oenv);
1214   }
1215   fprintf(stdout,"\n%-24s %10s\n",
1216           "Energy","dF = -kT ln < exp(-(EB-EA)/kT) >A");
1217   sum=0;
1218   beta = 1.0/(BOLTZ*reftemp);
1219   for(i=0; i<nset; i++) {
1220     if (gmx_strcasecmp(leg[i],enm[set[i]].name)!=0)
1221       fprintf(stderr,"\nWARNING energy set name mismatch %s!=%s\n",
1222               leg[i],enm[set[i]].name);
1223     for(j=0; j<nenergy; j++) {
1224       dE = eneset2[i][j] - edat->s[i].ener[j];
1225       sum += exp(-dE*beta);
1226       if (fp)
1227         fprintf(fp,"%10g %10g %10g\n",
1228                 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1229     }
1230     aver = -BOLTZ*reftemp*log(sum/nenergy);
1231     fprintf(stdout,"%-24s %10g\n",leg[i],aver);
1232   }
1233   if(fp) ffclose(fp);
1234   sfree(fr);
1235 }
1236
1237 int gmx_energy(int argc,char *argv[])
1238 {
1239   const char *desc[] = {
1240
1241     "g_energy extracts energy components or distance restraint",
1242     "data from an energy file. The user is prompted to interactively",
1243     "select the energy terms she wants.[PAR]",
1244
1245     "Average, RMSD and drift are calculated with full precision from the",
1246     "simulation (see printed manual). Drift is calculated by performing",
1247     "a LSQ fit of the data to a straight line. The reported total drift",
1248     "is the difference of the fit at the first and last point.",
1249     "An error estimate of the average is given based on a block averages",
1250     "over 5 blocks using the full precision averages. The error estimate",
1251     "can be performed over multiple block lengths with the options",
1252     "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1253     "Note that in most cases the energy files contains averages over all",
1254     "MD steps, or over many more points than the number of frames in",
1255     "energy file. This makes the g_energy statistics output more accurate",
1256     "than the xvg output. When exact averages are not present in the energy",
1257     "file the statistics mentioned above is simply over the single, per-frame",
1258     "energy values.[PAR]",
1259
1260     "The term fluctuation gives the RMSD around the LSQ fit.[PAR]",
1261
1262     "Some fluctuation-dependent properties can be calculated provided",
1263     "the correct energy terms are selected. The following properties",
1264     "will be computed:[BR]",
1265     "Property                        Energy terms needed[BR]",
1266     "---------------------------------------------------[BR]",
1267     "Heat capacity Cp (NPT sims):    Enthalpy, Temp     [BR]",
1268     "Heat capacity Cv (NVT sims):    Etot, Temp         [BR]",
1269     "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1270     "Isothermal compressibility:     Vol, Temp          [BR]",
1271     "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1272     "---------------------------------------------------[BR]",
1273     "You always need to set the number of molecules [TT]-nmol[tt], and,",
1274     "if you used constraints in your simulations you will need to give",
1275     "the number of constraints per molecule [TT]-nconstr[tt] in order to",
1276     "correct for this: (nconstr/2) kB is subtracted from the heat",
1277     "capacity in this case. For instance in the case of rigid water",
1278     "you need to give the value 3 to this option.[PAR]",
1279
1280     "When the [TT]-viol[tt] option is set, the time averaged",
1281     "violations are plotted and the running time-averaged and",
1282     "instantaneous sum of violations are recalculated. Additionally",
1283     "running time-averaged and instantaneous distances between",
1284     "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1285
1286     "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1287     "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1288     "The first two options plot the orientation, the last three the",
1289     "deviations of the orientations from the experimental values.",
1290     "The options that end on an 'a' plot the average over time",
1291     "as a function of restraint. The options that end on a 't'",
1292     "prompt the user for restraint label numbers and plot the data",
1293     "as a function of time. Option [TT]-odr[tt] plots the RMS",
1294     "deviation as a function of restraint.",
1295     "When the run used time or ensemble averaged orientation restraints,",
1296     "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1297     "not ensemble-averaged orientations and deviations instead of",
1298     "the time and ensemble averages.[PAR]",
1299
1300     "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1301     "tensor for each orientation restraint experiment. With option",
1302     "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1303
1304     "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1305     "difference with an ideal gas state: [BR]",
1306     "  Delta A = A(N,V,T) - A_idgas(N,V,T) = kT ln < e^(Upot/kT) >[BR]",
1307     "  Delta G = G(N,p,T) - G_idgas(N,p,T) = kT ln < e^(Upot/kT) >[BR]",
1308     "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1309     "the average is over the ensemble (or time in a trajectory).",
1310     "Note that this is in principle",
1311     "only correct when averaging over the whole (Boltzmann) ensemble",
1312     "and using the potential energy. This also allows for an entropy",
1313     "estimate using:[BR]",
1314     "  Delta S(N,V,T) = S(N,V,T) - S_idgas(N,V,T) = (<Upot> - Delta A)/T[BR]",
1315     "  Delta S(N,p,T) = S(N,p,T) - S_idgas(N,p,T) = (<Upot> + pV - Delta G)/T",
1316     "[PAR]",
1317
1318     "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1319     "difference is calculated dF = -kT ln < e ^ -(EB-EA)/kT >A ,",
1320     "where EA and EB are the energies from the first and second energy",
1321     "files, and the average is over the ensemble A. [BB]NOTE[bb] that",
1322     "the energies must both be calculated from the same trajectory."
1323
1324   };
1325   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE;
1326   static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
1327   static int  skip=0,nmol=1,nconstr=0,nbmin=5,nbmax=5;
1328   static real reftemp=300.0,ezero=0;
1329   t_pargs pa[] = {
1330     { "-fee",   FALSE, etBOOL,  {&bFee},
1331       "Do a free energy estimate" },
1332     { "-fetemp", FALSE, etREAL,{&reftemp},
1333       "Reference temperature for free energy calculation" },
1334     { "-zero", FALSE, etREAL, {&ezero},
1335       "Subtract a zero-point energy" },
1336     { "-sum",  FALSE, etBOOL, {&bSum},
1337       "Sum the energy terms selected rather than display them all" },
1338     { "-dp",   FALSE, etBOOL, {&bDp},
1339       "Print energies in high precision" },
1340     { "-nbmin", FALSE, etINT, {&nbmin},
1341       "Minimum number of blocks for error estimate" },
1342     { "-nbmax", FALSE, etINT, {&nbmax},
1343       "Maximum number of blocks for error estimate" },
1344     { "-mutot",FALSE, etBOOL, {&bMutot},
1345       "Compute the total dipole moment from the components" },
1346     { "-skip", FALSE, etINT,  {&skip},
1347       "Skip number of frames between data points" },
1348     { "-aver", FALSE, etBOOL, {&bPrAll},
1349       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1350     { "-nmol", FALSE, etINT,  {&nmol},
1351       "Number of molecules in your sample: the energies are divided by this number" },
1352     { "-nconstr",  FALSE, etINT,  {&nconstr},
1353       "Number of constraints per molecule. Necessary for calculating the heat capacity" },
1354     { "-fluc", FALSE, etBOOL, {&bFluct},
1355       "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1356     { "-orinst", FALSE, etBOOL, {&bOrinst},
1357       "Analyse instantaneous orientation data" },
1358     { "-ovec", FALSE, etBOOL, {&bOvec},
1359       "Also plot the eigenvectors with -oten" }
1360   };
1361   const char* drleg[] = {
1362     "Running average",
1363     "Instantaneous"
1364   };
1365   static const char *setnm[] = {
1366     "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1367     "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1368     "Volume",  "Pressure"
1369   };
1370
1371   FILE       *out,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1372   FILE       **drout;
1373   ener_file_t fp;
1374   int        timecheck=0;
1375   gmx_mtop_t mtop;
1376   gmx_localtop_t *top=NULL;
1377   t_inputrec ir;
1378   t_energy   **ee;
1379   enerdata_t edat;
1380   gmx_enxnm_t *enm=NULL;
1381   t_enxframe *frame,*fr=NULL;
1382   int        cur=0;
1383 #define NEXT (1-cur)
1384   int        nre,teller,teller_disre,nfr;
1385   gmx_large_int_t start_step;
1386   int        nor=0,nex=0,norfr=0,enx_i=0;
1387   real       start_t;
1388   real       *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1389   int        *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1390   int        nbounds=0,npairs;
1391   gmx_bool       bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN;
1392   gmx_bool       bFoundStart,bCont,bEDR,bVisco;
1393   double     sum,sumaver,sumt,ener,dbl;
1394   double     *time=NULL;
1395   real       Vaver;
1396   int        *set=NULL,i,j,k,nset,sss;
1397   gmx_bool       *bIsEner=NULL;
1398   char       **pairleg,**odtleg,**otenleg;
1399   char       **leg=NULL;
1400   char       **nms;
1401   char       *anm_j,*anm_k,*resnm_j,*resnm_k;
1402   int        resnr_j,resnr_k;
1403   const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1404   char       buf[256];
1405   output_env_t oenv;
1406   t_enxblock *blk=NULL;
1407   t_enxblock *blk_disre=NULL;
1408   int        ndisre=0;
1409
1410   t_filenm   fnm[] = {
1411     { efEDR, "-f",    NULL,      ffREAD  },
1412     { efEDR, "-f2",   NULL,      ffOPTRD },
1413     { efTPX, "-s",    NULL,      ffOPTRD },
1414     { efXVG, "-o",    "energy",  ffWRITE },
1415     { efXVG, "-viol", "violaver",ffOPTWR },
1416     { efXVG, "-pairs","pairs",   ffOPTWR },
1417     { efXVG, "-ora",  "orienta", ffOPTWR },
1418     { efXVG, "-ort",  "orientt", ffOPTWR },
1419     { efXVG, "-oda",  "orideva", ffOPTWR },
1420     { efXVG, "-odr",  "oridevr", ffOPTWR },
1421     { efXVG, "-odt",  "oridevt", ffOPTWR },
1422     { efXVG, "-oten", "oriten",  ffOPTWR },
1423     { efXVG, "-corr", "enecorr", ffOPTWR },
1424     { efXVG, "-vis",  "visco",   ffOPTWR },
1425     { efXVG, "-ravg", "runavgdf",ffOPTWR }
1426   };
1427 #define NFILE asize(fnm)
1428   int     npargs;
1429   t_pargs *ppa;
1430
1431   CopyRight(stderr,argv[0]);
1432   npargs = asize(pa);
1433   ppa    = add_acf_pargs(&npargs,pa);
1434   parse_common_args(&argc,argv,
1435                     PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1436                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1437
1438   bDRAll = opt2bSet("-pairs",NFILE,fnm);
1439   bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1440   bORA   = opt2bSet("-ora",NFILE,fnm);
1441   bORT   = opt2bSet("-ort",NFILE,fnm);
1442   bODA   = opt2bSet("-oda",NFILE,fnm);
1443   bODR   = opt2bSet("-odr",NFILE,fnm);
1444   bODT   = opt2bSet("-odt",NFILE,fnm);
1445   bORIRE = bORA || bORT || bODA || bODR || bODT;
1446   bOTEN  = opt2bSet("-oten",NFILE,fnm);
1447
1448   nset = 0;
1449
1450   snew(frame,2);
1451   fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1452   do_enxnms(fp,&nre,&enm);
1453
1454   Vaver = -1;
1455
1456   bVisco = opt2bSet("-vis",NFILE,fnm);
1457
1458   if (!bDisRe) {
1459     if (bVisco) {
1460       nset=asize(setnm);
1461       snew(set,nset);
1462       /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1463       for(j=0; j<nset; j++) {
1464         for(i=0; i<nre; i++) {
1465           if (strstr(enm[i].name,setnm[j])) {
1466             set[j]=i;
1467             break;
1468           }
1469         }
1470         if (i == nre) {
1471           if (gmx_strcasecmp(setnm[j],"Volume")==0) {
1472             printf("Enter the box volume (" unit_volume "): ");
1473             if(1 != scanf("%lf",&dbl))
1474             {
1475               gmx_fatal(FARGS,"Error reading user input");
1476             }
1477             Vaver = dbl;
1478           } else
1479             gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1480                         setnm[j]);
1481         }
1482       }
1483     }
1484     else {
1485       set=select_by_name(nre,enm,&nset);
1486     }
1487     /* Print all the different units once */
1488     sprintf(buf,"(%s)",enm[set[0]].unit);
1489     for(i=1; i<nset; i++) {
1490       for(j=0; j<i; j++) {
1491         if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1492           break;
1493         }
1494       }
1495       if (j == i) {
1496         strcat(buf,", (");
1497         strcat(buf,enm[set[i]].unit);
1498         strcat(buf,")");
1499       }
1500     }
1501     out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1502                  oenv);
1503
1504     snew(leg,nset+1);
1505     for(i=0; (i<nset); i++)
1506       leg[i] = enm[set[i]].name;
1507     if (bSum) {
1508       leg[nset]=strdup("Sum");
1509       xvgr_legend(out,nset+1,(const char**)leg,oenv);
1510     }
1511     else
1512       xvgr_legend(out,nset,(const char**)leg,oenv);
1513
1514     snew(bIsEner,nset);
1515     for(i=0; (i<nset); i++) {
1516       bIsEner[i] = FALSE;
1517       for (j=0; (j <= F_ETOT); j++)
1518         bIsEner[i] = bIsEner[i] ||
1519           (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1520     }
1521
1522     if (bPrAll && nset > 1) {
1523         gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1524     }
1525
1526     time = NULL;
1527
1528     if (bORIRE || bOTEN)
1529       get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1530
1531     if (bORIRE) {
1532       if (bOrinst)
1533         enx_i = enxORI;
1534       else
1535         enx_i = enxOR;
1536
1537       if (bORA || bODA)
1538         snew(orient,nor);
1539       if (bODR)
1540         snew(odrms,nor);
1541       if (bORT || bODT) {
1542         fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1543         fprintf(stderr,"End your selection with 0\n");
1544         j = -1;
1545         orsel = NULL;
1546         do {
1547           j++;
1548           srenew(orsel,j+1);
1549           if(1 != scanf("%d",&(orsel[j])))
1550           {
1551             gmx_fatal(FARGS,"Error reading user input");
1552           }
1553         } while (orsel[j] > 0);
1554         if (orsel[0] == -1) {
1555           fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1556           norsel = nor;
1557           srenew(orsel,nor);
1558           for(i=0; i<nor; i++)
1559             orsel[i] = i;
1560         } else {
1561           /* Build the selection */
1562           norsel=0;
1563           for(i=0; i<j; i++) {
1564             for(k=0; k<nor; k++)
1565               if (or_label[k] == orsel[i]) {
1566                 orsel[norsel] = k;
1567                 norsel++;
1568                 break;
1569               }
1570             if (k == nor)
1571               fprintf(stderr,"Orientation restraint label %d not found\n",
1572                       orsel[i]);
1573           }
1574         }
1575         snew(odtleg,norsel);
1576         for(i=0; i<norsel; i++) {
1577           snew(odtleg[i],256);
1578           sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1579         }
1580         if (bORT) {
1581           fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1582                         "Time (ps)","",oenv);
1583           if (bOrinst)
1584             fprintf(fort,"%s",orinst_sub);
1585           xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
1586         }
1587         if (bODT) {
1588           fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1589                         "Orientation restraint deviation",
1590                         "Time (ps)","",oenv);
1591           if (bOrinst)
1592             fprintf(fodt,"%s",orinst_sub);
1593           xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
1594         }
1595       }
1596     }
1597     if (bOTEN) {
1598       foten=xvgropen(opt2fn("-oten",NFILE,fnm),
1599                      "Order tensor","Time (ps)","",oenv);
1600       snew(otenleg,bOvec ? nex*12 : nex*3);
1601       for(i=0; i<nex; i++) {
1602         for(j=0; j<3; j++) {
1603           sprintf(buf,"eig%d",j+1);
1604           otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
1605         }
1606         if (bOvec) {
1607           for(j=0; j<9; j++) {
1608             sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
1609             otenleg[12*i+3+j] = strdup(buf);
1610           }
1611         }
1612       }
1613       xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
1614     }
1615   }
1616   else {
1617     nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
1618                        &mtop,&top,&ir);
1619     snew(violaver,npairs);
1620     out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
1621                  "Time (ps)","nm",oenv);
1622     xvgr_legend(out,2,drleg,oenv);
1623     if (bDRAll) {
1624       fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
1625                         "Time (ps)","Distance (nm)",oenv);
1626       if (output_env_get_print_xvgr_codes(oenv))
1627         fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
1628                 ir.dr_tau);
1629     }
1630   }
1631
1632   /* Initiate energies and set them to zero */
1633   edat.nsteps  = 0;
1634   edat.npoints = 0;
1635   edat.nframes = 0;
1636   edat.step    = NULL;
1637   edat.steps   = NULL;
1638   edat.points  = NULL;
1639   snew(edat.s,nset);
1640
1641   /* Initiate counters */
1642   teller       = 0;
1643   teller_disre = 0;
1644   bFoundStart  = FALSE;
1645   start_step   = 0;
1646   start_t      = 0;
1647   do {
1648     /* This loop searches for the first frame (when -b option is given),
1649      * or when this has been found it reads just one energy frame
1650      */
1651     do {
1652       bCont = do_enx(fp,&(frame[NEXT]));
1653
1654       if (bCont) {
1655         timecheck = check_times(frame[NEXT].t);
1656       }
1657     } while (bCont && (timecheck < 0));
1658
1659     if ((timecheck == 0) && bCont) {
1660       /* We read a valid frame, so we can use it */
1661       fr = &(frame[NEXT]);
1662
1663       if (fr->nre > 0) {
1664         /* The frame contains energies, so update cur */
1665         cur  = NEXT;
1666
1667                 if (edat.nframes % 1000 == 0)
1668             {
1669                 srenew(edat.step,edat.nframes+1000);
1670                 srenew(edat.steps,edat.nframes+1000);
1671                 srenew(edat.points,edat.nframes+1000);
1672                 for(i=0; i<nset; i++)
1673                 {
1674                     srenew(edat.s[i].ener,edat.nframes+1000);
1675                     srenew(edat.s[i].es  ,edat.nframes+1000);
1676                 }
1677             }
1678
1679                 nfr = edat.nframes;
1680             edat.step[nfr] = fr->step;
1681
1682             if (!bFoundStart)
1683             {
1684                 bFoundStart = TRUE;
1685                 /* Initiate the previous step data */
1686                 start_step = fr->step;
1687                 start_t    = fr->t;
1688                 /* Initiate the energy sums */
1689                 edat.steps[nfr]  = 1;
1690                 edat.points[nfr] = 1;
1691                 for(i=0; i<nset; i++)
1692                 {
1693                     sss = set[i];
1694                     edat.s[i].es[nfr].sum  = fr->ener[sss].e;
1695                     edat.s[i].es[nfr].sum2 = 0;
1696                 }
1697                 edat.nsteps  = 1;
1698                 edat.npoints = 1;
1699             }
1700             else
1701             {
1702                 edat.steps[nfr] = fr->nsteps;
1703                 {
1704                     if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
1705                     {
1706                         if (fr->nsum <= 1)
1707                         {
1708                             edat.points[nfr] = 1;
1709                             for(i=0; i<nset; i++)
1710                             {
1711                                 sss = set[i];
1712                                 edat.s[i].es[nfr].sum  = fr->ener[sss].e;
1713                                 edat.s[i].es[nfr].sum2 = 0;
1714                             }
1715                             edat.npoints += 1;
1716                         }
1717                         else
1718                         {
1719                             edat.points[nfr] = fr->nsum;
1720                             for(i=0; i<nset; i++)
1721                             {
1722                                 sss = set[i];
1723                                 edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
1724                                 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
1725                             }
1726                             edat.npoints += fr->nsum;
1727                         }
1728                     }
1729                     else
1730                     {
1731                         /* The interval does not match fr->nsteps:
1732                          * can not do exact averages.
1733                          */
1734                         edat.npoints = 0;
1735                     }
1736                     edat.nsteps = fr->step - start_step + 1;
1737                 }
1738             }
1739             for(i=0; i<nset; i++)
1740             {
1741                 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
1742             }
1743       }
1744       /*
1745        * Define distance restraint legends. Can only be done after
1746        * the first frame has been read... (Then we know how many there are)
1747        */
1748       blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
1749       if (bDisRe && bDRAll && !leg && blk_disre)
1750       {
1751           t_iatom   *fa;
1752           t_iparams *ip;
1753
1754           fa = top->idef.il[F_DISRES].iatoms;
1755           ip = top->idef.iparams;
1756           if (blk_disre->nsub != 2 ||
1757               (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
1758           {
1759               gmx_incons("Number of disre sub-blocks not equal to 2");
1760           }
1761
1762           ndisre=blk_disre->sub[0].nr ;
1763           if (ndisre != top->idef.il[F_DISRES].nr/3)
1764           {
1765               gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
1766                         ndisre,top->idef.il[F_DISRES].nr/3);
1767           }
1768           snew(pairleg,ndisre);
1769           for(i=0; i<ndisre; i++)
1770           {
1771               snew(pairleg[i],30);
1772               j=fa[3*i+1];
1773               k=fa[3*i+2];
1774               gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
1775               gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
1776               sprintf(pairleg[i],"%d %s %d %s (%d)",
1777                       resnr_j,anm_j,resnr_k,anm_k,
1778                       ip[fa[3*i]].disres.label);
1779           }
1780           set=select_it(ndisre,pairleg,&nset);
1781           snew(leg,2*nset);
1782           for(i=0; (i<nset); i++)
1783           {
1784               snew(leg[2*i],32);
1785               sprintf(leg[2*i],  "a %s",pairleg[set[i]]);
1786               snew(leg[2*i+1],32);
1787               sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
1788           }
1789           xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
1790       }
1791
1792       /*
1793        * Store energies for analysis afterwards...
1794        */
1795       if (!bDisRe && (fr->nre > 0)) {
1796         if (edat.nframes % 1000 == 0) {
1797           srenew(time,edat.nframes+1000);
1798         }
1799         time[edat.nframes] = fr->t;
1800         edat.nframes++;
1801       }
1802       /*
1803        * Printing time, only when we do not want to skip frames
1804        */
1805       if (!skip || teller % skip == 0) {
1806         if (bDisRe) {
1807           /*******************************************
1808            * D I S T A N C E   R E S T R A I N T S
1809            *******************************************/
1810           if (ndisre > 0)
1811           {
1812 #ifndef GMX_DOUBLE
1813             float *disre_rt =     blk_disre->sub[0].fval;
1814             float *disre_rm3tav = blk_disre->sub[1].fval;
1815 #else
1816             double *disre_rt =     blk_disre->sub[0].dval;
1817             double *disre_rm3tav = blk_disre->sub[1].dval;
1818 #endif
1819
1820             print_time(out,fr->t);
1821             if (violaver == NULL)
1822               snew(violaver,ndisre);
1823
1824             /* Subtract bounds from distances, to calculate violations */
1825             calc_violations(disre_rt, disre_rm3tav,
1826                             nbounds,pair,bounds,violaver,&sumt,&sumaver);
1827
1828             fprintf(out,"  %8.4f  %8.4f\n",sumaver,sumt);
1829             if (bDRAll) {
1830               print_time(fp_pairs,fr->t);
1831               for(i=0; (i<nset); i++) {
1832                 sss=set[i];
1833                 fprintf(fp_pairs,"  %8.4f", mypow(disre_rm3tav[sss],minthird));
1834                 fprintf(fp_pairs,"  %8.4f", disre_rt[sss]);
1835               }
1836               fprintf(fp_pairs,"\n");
1837             }
1838             teller_disre++;
1839           }
1840         }
1841         /*******************************************
1842          * E N E R G I E S
1843          *******************************************/
1844         else {
1845           if (fr->nre > 0) {
1846             if (bPrAll)
1847             {
1848                 /* We skip frames with single points (usually only the first frame),
1849                  * since they would result in an average plot with outliers.
1850                  */
1851                 if (fr->nsum > 1) {
1852                     print_time(out,fr->t);
1853                      print1(out,bDp,fr->ener[set[0]].e);
1854                      print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
1855                      print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
1856                      fprintf(out,"\n");
1857                 }
1858             }
1859             else
1860             {
1861                 print_time(out,fr->t);
1862                 if (bSum)
1863                 {
1864                     sum = 0;
1865                     for(i=0; i<nset; i++)
1866                     {
1867                         sum += fr->ener[set[i]].e;
1868                     }
1869                     print1(out,bDp,sum/nmol-ezero);
1870                 }
1871                 else
1872                 {
1873                     for(i=0; (i<nset); i++)
1874                     {
1875                         if (bIsEner[i])
1876                         {
1877                             print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
1878                         }
1879                         else
1880                         {
1881                             print1(out,bDp,fr->ener[set[i]].e);
1882                         }
1883                     }
1884                 }
1885                 fprintf(out,"\n");
1886             }
1887           }
1888 #if 0
1889           /* we first count the blocks that have id 0: the orire blocks */
1890           block_orire=0;
1891           for(b=0;b<fr->nblock;b++)
1892           {
1893               if (fr->block[b].id == mde_block_type_orire)
1894                   nblock_orire++;
1895           }
1896 #endif
1897           blk = find_block_id_enxframe(fr, enx_i, NULL);
1898           if (bORIRE && blk)
1899           {
1900 #ifndef GMX_DOUBLE
1901               xdr_datatype dt=xdr_datatype_float;
1902 #else
1903               xdr_datatype dt=xdr_datatype_double;
1904 #endif
1905               real *vals;
1906
1907               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
1908               {
1909                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
1910               }
1911 #ifndef GMX_DOUBLE
1912               vals=blk->sub[0].fval;
1913 #else
1914               vals=blk->sub[0].dval;
1915 #endif
1916
1917               if (blk->sub[0].nr != (size_t)nor)
1918                   gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
1919               if (bORA || bODA)
1920               {
1921                   for(i=0; i<nor; i++)
1922                       orient[i] += vals[i];
1923               }
1924               if (bODR)
1925               {
1926                   for(i=0; i<nor; i++)
1927                       odrms[i] += sqr(vals[i]-oobs[i]);
1928               }
1929               if (bORT)
1930               {
1931                   fprintf(fort,"  %10f",fr->t);
1932                   for(i=0; i<norsel; i++)
1933                       fprintf(fort," %g",vals[orsel[i]]);
1934                   fprintf(fort,"\n");
1935               }
1936               if (bODT)
1937               {
1938                   fprintf(fodt,"  %10f",fr->t);
1939                   for(i=0; i<norsel; i++)
1940                       fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
1941                   fprintf(fodt,"\n");
1942               }
1943               norfr++;
1944           }
1945           blk = find_block_id_enxframe(fr, enxORT, NULL);
1946           if (bOTEN && blk)
1947           {
1948 #ifndef GMX_DOUBLE
1949               xdr_datatype dt=xdr_datatype_float;
1950 #else
1951               xdr_datatype dt=xdr_datatype_double;
1952 #endif
1953               real *vals;
1954
1955               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
1956                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
1957 #ifndef GMX_DOUBLE
1958               vals=blk->sub[0].fval;
1959 #else
1960               vals=blk->sub[0].dval;
1961 #endif
1962
1963               if (blk->sub[0].nr != (size_t)(nex*12))
1964                   gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
1965                             blk->sub[0].nr/12, nex);
1966               fprintf(foten,"  %10f",fr->t);
1967               for(i=0; i<nex; i++)
1968                   for(j=0; j<(bOvec?12:3); j++)
1969                       fprintf(foten," %g",vals[i*12+j]);
1970               fprintf(foten,"\n");
1971           }
1972         }
1973       }
1974     }
1975   } while (bCont && (timecheck == 0));
1976
1977   fprintf(stderr,"\n");
1978   close_enx(fp);
1979
1980   ffclose(out);
1981
1982   if (bDRAll)
1983     ffclose(fp_pairs);
1984
1985   if (bORT)
1986     ffclose(fort);
1987   if (bODT)
1988     ffclose(fodt);
1989   if (bORA) {
1990     out = xvgropen(opt2fn("-ora",NFILE,fnm),
1991                    "Average calculated orientations",
1992                    "Restraint label","",oenv);
1993     if (bOrinst)
1994       fprintf(out,"%s",orinst_sub);
1995     for(i=0; i<nor; i++)
1996       fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr);
1997     ffclose(out);
1998   }
1999   if (bODA) {
2000     out = xvgropen(opt2fn("-oda",NFILE,fnm),
2001                    "Average restraint deviation",
2002                    "Restraint label","",oenv);
2003     if (bOrinst)
2004       fprintf(out,"%s",orinst_sub);
2005     for(i=0; i<nor; i++)
2006       fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr-oobs[i]);
2007     ffclose(out);
2008   }
2009   if (bODR) {
2010     out = xvgropen(opt2fn("-odr",NFILE,fnm),
2011                    "RMS orientation restraint deviations",
2012                    "Restraint label","",oenv);
2013     if (bOrinst)
2014       fprintf(out,"%s",orinst_sub);
2015     for(i=0; i<nor; i++)
2016       fprintf(out,"%5d  %g\n",or_label[i],sqrt(odrms[i]/norfr));
2017     ffclose(out);
2018   }
2019   if (bOTEN)
2020     ffclose(foten);
2021
2022   if (bDisRe) {
2023     analyse_disre(opt2fn("-viol",NFILE,fnm),
2024                   teller_disre,violaver,bounds,index,pair,nbounds,oenv);
2025   } else {
2026     analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
2027                  bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
2028                  bVisco,opt2fn("-vis",NFILE,fnm),
2029                  nmol,nconstr,start_step,start_t,frame[cur].step,frame[cur].t,
2030                  time,reftemp,&edat,
2031                  nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
2032                  oenv);
2033   }
2034   if (opt2bSet("-f2",NFILE,fnm)) {
2035     fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
2036         reftemp, nset, set, leg, &edat, time ,oenv);
2037   }
2038
2039   {
2040     const char *nxy = "-nxy";
2041
2042     do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
2043     do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
2044     do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
2045     do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
2046     do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
2047     do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
2048     do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
2049     do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
2050   }
2051   thanx(stderr);
2052
2053   return 0;
2054 }