src/tools/g_energy.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_g_energy_c = "$Id$";
  30
  31 #include <string.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "fatal.h"
  35 #include "vec.h"
  36 #include "smalloc.h"
  37 #include "enxio.h"
  38 #include "statutil.h"
  39 #include "assert.h"
  40 #include "names.h"
  41 #include "copyrite.h"
  42 #include "macros.h"
  43 #include "xvgr.h"
  44 #include "gstat.h"
  45 #include "physics.h"
  46 #include "tpxio.h"
  47
  48 static real       minthird=-1.0/3.0,minsixth=-1.0/6.0;
  49
  50 double mypow(double x,double y)
  51 {
  52   if (x > 0)
  53     return pow(x,y);
  54   else
  55     return 0.0;
  56 }
  57
  58 static int *select_it(int nre,char *nm[],int *nset)
  59 {
  60   bool *bE;
  61   int  n,k,j,i;
  62   int  *set;
  63   bool bVerbose = TRUE;
  64
  65   if ((getenv("VERBOSE")) != NULL)
  66     bVerbose = FALSE;
  67
  68   printf("Select the terms you want from the following list\n");
  69   printf("End your selection with 0\n");
  70
  71   if ( bVerbose ) {
  72     for(k=0; (k<nre); ) {
  73       for(j=0; (j<4) && (k<nre); j++,k++)
  74         printf(" %3d=%14s",k+1,nm[k]);
  75       printf("\n");
  76     }
  77   }
  78
  79   snew(bE,nre);
  80   do {
  81     scanf("%d",&n);
  82     if ((n>0) && (n<=nre))
  83       bE[n-1]=TRUE;
  84   } while (n != 0);
  85
  86   snew(set,nre);
  87   for(i=(*nset)=0; (i<nre); i++)
  88     if (bE[i])
  89       set[(*nset)++]=i;
  90
  91   sfree(bE);
  92
  93   return set;
  94 }
  95
  96 int get_bounds(char *topnm,real **bounds,int **index,int **dr_pair,int *npairs,
  97                t_topology *top,t_inputrec *ir)
  98 {
  99   t_functype *functype;
 100   t_iparams  *ip;
 101   int        natoms,i,j,k,type,ftype,natom;
 102   t_ilist    *disres;
 103   t_iatom    *iatom;
 104   real       *b,t;
 105   int        *ind,*pair;
 106   int        nb,index1;
 107   matrix     box;
 108
 109   read_tpx(topnm,&i,&t,&t,ir,box,&natoms,NULL,NULL,NULL,top);
 110
 111   functype = top->idef.functype;
 112   ip       = top->idef.iparams;
 113
 114   /* Count how many distance restraint there are... */
 115   nb=top->idef.il[F_DISRES].nr;
 116   if (nb == 0)
 117     fatal_error(0,"No distance restraints in topology!\n");
 118
 119   /* Allocate memory */
 120   snew(b,nb);
 121   snew(ind,nb);
 122   snew(pair,nb+1);
 123
 124   /* Fill the bound array */
 125   nb=0;
 126   for(i=0; (i<top->idef.ntypes); i++) {
 127     ftype = functype[i];
 128     if (ftype == F_DISRES) {
 129       index1 = ip[i].disres.index;
 130       b[nb]   = ip[i].disres.up1;
 131       ind[nb] = index1;
 132       nb++;
 133     }
 134   }
 135   *bounds = b;
 136
 137   /* Fill the index array */
 138   index1  = -1;
 139   disres  = &(top->idef.il[F_DISRES]);
 140   iatom   = disres->iatoms;
 141   for(i=j=k=0; (i<disres->nr); ) {
 142     type  = iatom[i];
 143     ftype = top->idef.functype[type];
 144     natom = interaction_function[ftype].nratoms+1;
 145     if (index1 != top->idef.iparams[type].disres.index) {
 146       pair[j] = k;
 147       index1  = top->idef.iparams[type].disres.index;
 148       j ++;
 149     }
 150     k++;
 151     i += natom;
 152   }
 153   pair[j]  = k;
 154   *npairs = k;
 155   assert(j == nb);
 156
 157   *index   = ind;
 158   *dr_pair = pair;
 159
 160   return nb;
 161 }
 162
 163 void calc_violations(real rt[],real rav3[],int nb,int index[],
 164                      real bounds[],real *viol,double *st,double *sa)
 165 {
 166   const   real sixth=1.0/6.0;
 167   int     i,j;
 168   double  rsum,rav,sumaver,sumt;
 169
 170   sumaver = 0;
 171   sumt    = 0;
 172   for(i=0; (i<nb); i++) {
 173     rsum = 0.0;
 174     rav  = 0.0;
 175     for(j=index[i]; (j<index[i+1]); j++) {
 176       if (viol)
 177         viol[j] += mypow(rt[j],-3.0);
 178       rav     += sqr(rav3[j]);
 179       rsum    += mypow(rt[j],-6);
 180     }
 181     rsum    = max(0.0,mypow(rsum,-sixth)-bounds[i]);
 182     rav     = max(0.0,mypow(rav, -sixth)-bounds[i]);
 183
 184     sumt    += rsum;
 185     sumaver += rav;
 186   }
 187   *st = sumt;
 188   *sa = sumaver;
 189 }
 190
 191 static void analyse_disre(char *voutfn,    int nframes,
 192                           real violaver[], real bounds[], int index[],
 193                           int pair[],      int nbounds)
 194 {
 195   FILE   *vout;
 196   double sum,sumt,sumaver;
 197   int    i,j;
 198
 199   /* Subtract bounds from distances, to calculate violations */
 200   calc_violations(violaver,violaver,
 201                   nbounds,pair,bounds,NULL,&sumt,&sumaver);
 202
 203 #ifdef DEBUG
 204   fprintf(stderr,"\nSum of violations averaged over simulation: %g nm\n",
 205           sumaver);
 206   fprintf(stderr,"Largest violation averaged over simulation: %g nm\n\n",
 207           sumt);
 208 #endif
 209   vout=xvgropen(voutfn,"r\\S-3\\N average violations","DR Index","nm");
 210   sum  = 0.0;
 211   sumt = 0.0;
 212   for(i=0; (i<nbounds); i++) {
 213     /* Do ensemble averaging */
 214     sumaver = 0;
 215     for(j=pair[i]; (j<pair[i+1]); j++)
 216       sumaver += sqr(violaver[j]/nframes);
 217     sumaver = max(0.0,mypow(sumaver,minsixth)-bounds[i]);
 218
 219     sumt   += sumaver;
 220     sum     = max(sum,sumaver);
 221     fprintf(vout,"%10d  %10.5e\n",index[i],sumaver);
 222   }
 223 #ifdef DEBUG
 224   for(j=0; (j<dr.ndr); j++)
 225     fprintf(vout,"%10d  %10.5e\n",j,mypow(violaver[j]/nframes,minthird));
 226 #endif
 227   ffclose(vout);
 228
 229   fprintf(stderr,"\nSum of violations averaged over simulation: %g nm\n",sumt);
 230   fprintf(stderr,"Largest violation averaged over simulation: %g nm\n\n",sum);
 231
 232   xvgr_file(voutfn,"-graphtype bar");
 233 }
 234
 235 static void einstein_visco(char *fn,int nsets,int nframes,real **set,
 236                            real V,real T,real time[])
 237 {
 238   FILE *fp;
 239   real fac,**integral,dt,di,r_nf2;
 240   int  i,j,m,nf2;
 241
 242   if (nframes < 1)
 243     return;
 244
 245   dt  = (time[1]-time[0]);
 246
 247   /* Store the average in integral[nsets] */
 248   nf2 = nframes/2;
 249   snew(integral,nsets+1);
 250   for(m=0; (m<nsets+1); m++)
 251     snew(integral[m],nf2);
 252
 253   /* Use trapezium rule for integration and average over time origins */
 254   for(j=0; (j<nf2); j++) {
 255     for(i=1; (i<nf2); i++) {
 256       for(m=0; (m<nsets); m++) {
 257         di   = 0.5*dt*(set[m][i+j]+set[m][i-1+j]);
 258
 259         integral[m][i]     += di;
 260         integral[nsets][i] += di/nsets;
 261       }
 262     }
 263   }
 264
 265   fp=xvgropen(fn,"Shear viscosity using Einstein relation","Time (ps)","cp");
 266   fac   = (V*1e-26)/(2*BOLTZMANN*T);
 267   r_nf2 = 1.0/nf2;
 268   for(i=0; (i<nf2); i++) {
 269     fprintf(fp,"%10g",time[i]);
 270     for(m=0; (m<=nsets); m++)
 271       fprintf(fp,"  %10g",fac*sqr(r_nf2*integral[m][i]));
 272     fprintf(fp,"\n");
 273   }
 274   fclose(fp);
 275   for(m=0; (m<nsets+1); m++)
 276     sfree(integral[m]);
 277   sfree(integral);
 278 }
 279
 280 static void analyse_ener(bool bCorr,char *corrfn,
 281                          bool bGibbs,bool bSum,bool bFluct,
 282                          bool bVisco,char *visfn,int  nmol,int ndf,
 283                          int  oldstep,real oldt,int step,real t,
 284                          real time[], real reftemp,
 285                          t_energy oldee[],t_energy ee[],
 286                          int nset,int set[],int nenergy,real **eneset,
 287                          char *leg[])
 288 {
 289   FILE *fp;
 290   /* Check out the printed manual for equations! */
 291   real Dt,a,b,aver,avertot,stddev,delta_t,sigma,totaldrift;
 292   real xxx,integral,intBulk;
 293   real sfrac,oldfrac,diffsum,diffav,fstep,pr_aver,pr_stddev,fluct2;
 294   double beta=0,expE,expEtot,*gibbs=NULL;
 295   int  nsteps,iset;
 296   real x1m,x1mk,Temp=-1,Pres=-1,VarV=-1,Vaver=-1,VarT=-1;
 297   int  i,j,m,k;
 298   char buf[256];
 299
 300   nsteps  = step - oldstep;
 301 #ifdef DEBUG
 302   fprintf(stderr,"oldstep: %d, oldt: %g, step: %d, t: %g, nenergy: %d\n",
 303           oldstep,oldt,step,t,nenergy);
 304 #endif
 305   if (nsteps < 2) {
 306     fprintf(stderr,"Not enough steps (%d) for statistics\n",nsteps);
 307   }
 308   else {
 309     /* Calculate the time difference */
 310     delta_t = t - oldt;
 311
 312     fprintf(stderr,"\nStatistics over %d steps [ %.4f thru %.4f ps ], %d data sets\n\n",
 313             nsteps,oldt,t,nset);
 314
 315     fprintf(stderr,"%-24s %10s %10s %10s %10s %10s",
 316             "Energy","Average","RMSD","Fluct.","Drift","Tot-Drift");
 317     if (bGibbs)
 318       fprintf(stderr,"  %10s\n","-ln<e^(E/kT)>*kT");
 319     else
 320       fprintf(stderr,"\n");
 321     fprintf(stderr,"-------------------------------------------------------------------------------\n");
 322
 323     /* Initiate locals, only used with -sum */
 324     avertot=0;
 325     expEtot=0;
 326     if (bGibbs) {
 327       beta = 1.0/(BOLTZ*reftemp);
 328       snew(gibbs,nset);
 329     }
 330     for(i=0; (i<nset); i++) {
 331       iset = set[i];
 332       m      = oldstep;
 333       k      = nsteps;
 334 #ifdef DEBUG
 335       fprintf(stderr,"sum: %g, oldsum: %g, k: %d\n",
 336               ee[iset].esum,oldee[iset].esum,k);
 337 #endif
 338       aver   = (ee[iset].esum  - oldee[iset].esum)/k;
 339       fstep  = ((real) m) * ((real) (m+k))/((real) k);
 340       x1m    = (m > 0) ? oldee[iset].esum/m : 0.0;
 341       x1mk   = ee[iset].esum/(m+k);
 342       xxx    = sqr(x1m - x1mk);
 343       sigma  = ee[iset].eav - oldee[iset].eav - xxx * fstep;
 344       if((sigma/k)<GMX_REAL_EPS)
 345         sigma=0;
 346       stddev = sqrt(sigma/k);
 347
 348       if (bSum)
 349         avertot+=aver;
 350       if (bGibbs) {
 351         expE = 0;
 352         for(j=0; (j<nenergy); j++) {
 353           expE += exp(beta*(eneset[i][j]-aver)/nmol);
 354         }
 355         if (bSum)
 356           expEtot+=expE/nenergy;
 357
 358         gibbs[i] = log(expE/nenergy)/beta + aver/nmol;
 359       }
 360       if (strstr(leg[i],"empera") != NULL) {
 361         VarT = sqr(stddev);
 362         Temp = aver;
 363       } else if (strstr(leg[i],"olum") != NULL) {
 364         VarV = sqr(stddev);
 365         Vaver= aver;
 366       } else if (strstr(leg[i],"essure") != NULL) {
 367         Pres = aver;
 368       }
 369       if (iset < F_TEMP) {
 370         pr_aver   = aver/nmol;
 371         pr_stddev = stddev/nmol;
 372       }
 373       else {
 374         pr_aver   = aver;
 375         pr_stddev = stddev;
 376       }
 377       lsq_y_ax_b(nenergy,time,eneset[i],&a,&b);
 378       if(fabs(a)<GMX_REAL_EPS)
 379         a=0;
 380       totaldrift = a * delta_t * (nsteps+1)/nsteps;
 381       fluct2 = sqr(pr_stddev) - sqr(totaldrift)/12;
 382       if (fluct2 < 0)
 383         fluct2 = 0;
 384       fprintf(stderr,"%-24s %10g %10g %10g %10g %10g",
 385               leg[i],pr_aver,pr_stddev,sqrt(fluct2),a,totaldrift);
 386       if (bGibbs)
 387         fprintf(stderr,"  %10g\n",gibbs[i]);
 388       else
 389         fprintf(stderr,"\n");
 390       if (bFluct) {
 391         for(j=0; (j<nenergy); j++)
 392           eneset[i][j] -= aver;
 393       }
 394     }
 395     if (bSum) {
 396       fprintf(stderr,"%-24s %10g %10s %10s %10s %10s",
 397               "Total",avertot/nmol,"--","--","--","--");
 398       /* pr_aver,pr_stddev,a,totaldrift */
 399       if (bGibbs)
 400         fprintf(stderr,"  %10g  %10g\n",
 401                 log(expEtot)/beta + avertot/nmol,log(expEtot)/beta);
 402       else
 403         fprintf(stderr,"\n");
 404     }
 405     if (Temp != -1) {
 406       real factor;
 407
 408       factor = nmol*ndf*VarT/(3.0*sqr(Temp));
 409       fprintf(stderr,"Heat Capacity Cv:   %10g J/mol K (factor = %g)\n",
 410               1000*BOLTZ/(2.0/3.0 - factor),factor);
 411     }
 412     if ((VarV != -1) && (Temp != -1)) {
 413       real tmp = VarV/(Vaver*BOLTZ*Temp*PRESFAC);
 414
 415       fprintf(stderr,"Isothermal Compressibility: %10g /bar\n",tmp);
 416       fprintf(stderr,"Adiabatic bulk modulus:     %10g  bar\n",1.0/tmp);
 417     }
 418     /* Do correlation function */
 419     Dt = delta_t/nenergy;
 420     if (bVisco) {
 421       char *leg[] = { "Shear", "Bulk" };
 422       real factor;
 423
 424       /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
 425
 426       /* Symmetrise tensor! (and store in first three elements)
 427        * And subtract average pressure!
 428        */
 429       for(i=0; (i<nenergy); i++) {
 430         eneset[0][i] = 0.5*(eneset[1][i]+eneset[3][i]);
 431         eneset[1][i] = 0.5*(eneset[2][i]+eneset[6][i]);
 432         eneset[2][i] = 0.5*(eneset[5][i]+eneset[7][i]);
 433         eneset[11][i] -= Pres;
 434       }
 435
 436       einstein_visco("evisje.xvg",3,nenergy,eneset,Vaver,Temp,time);
 437
 438       /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
 439       /* Do it for shear viscosity */
 440       strcpy(buf,"Shear Viscosity");
 441       low_do_autocorr(corrfn,buf,nenergy,3,(nenergy+1)/2,eneset,Dt,
 442                       eacNormal,1,TRUE,TRUE,FALSE,FALSE,0.0,0.0,0);
 443
 444       /* Now for bulk viscosity */
 445       strcpy(buf,"Bulk Viscosity");
 446       low_do_autocorr(corrfn,buf,nenergy,1,(nenergy+1)/2,&(eneset[11]),Dt,
 447                       eacNormal,1,TRUE,TRUE,FALSE,FALSE,0.0,0.0,0);
 448
 449       factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
 450       fp=xvgropen(visfn,buf,"Time (ps)","\\8h\\4 (cp)");
 451       xvgr_legend(fp,asize(leg),leg);
 452
 453       /* Use trapezium rule for integration */
 454       integral = 0;
 455       intBulk  = 0;
 456       for(i=1; (i<nenergy/2); i++) {
 457         integral += 0.5*(eneset[0][i-1]  + eneset[0][i])*factor;
 458         intBulk  += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
 459         fprintf(fp,"%10g  %10g  %10g\n",(i*Dt),integral,intBulk);
 460       }
 461       fclose(fp);
 462     }
 463     else if (bCorr) {
 464       if (bFluct)
 465         strcpy(buf,"Autocorrelation of Energy Fluctuations");
 466       else
 467         strcpy(buf,"Energy Autocorrelation");
 468       do_autocorr(corrfn,buf,nenergy,nset,eneset,(delta_t/nenergy),
 469                   eacNormal,FALSE);
 470     }
 471   }
 472 }
 473
 474 void print_one(FILE *fp,bool bDp,real e)
 475 {
 476   if (bDp)
 477     fprintf(fp,"  %16g",e);
 478   else
 479     fprintf(fp,"  %10g",e);
 480
 481 }
 482
 483 int main(int argc,char *argv[])
 484 {
 485   static char *desc[] = {
 486
 487     "g_energy extracts energy components or distance restraint",
 488     "data from an energy file. The user is prompted to interactively",
 489     "select the energy terms she wants.[PAR]",
 490
 491     "When the [TT]-viol[tt] option is set, the time averaged",
 492     "violations are plotted and the running time-averaged and",
 493     "instantaneous sum of violations are recalculated. Additionally",
 494     "running time-averaged and instantaneous distances between",
 495     "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
 496
 497     "Average and RMSD are calculated with full precision from the",
 498     "simulation (see printed manual). Drift is calculated by performing",
 499     "a LSQ fit of the data to a straight line. Total drift is drift",
 500     "multiplied by total time.[PAR]",
 501
 502     "With [TT]-G[tt] a Gibbs free energy estimate is calculated using",
 503     "the formula: G = -ln < e ^ (E/kT) > * kT, where k is Boltzmann's",
 504     "constant, T is set by [TT]-Gtemp[tt] and the average is over the",
 505     "ensemble (or time in a trajectory). Note that this is in principle",
 506     "only correct when averaging over the whole (Boltzmann) ensemble",
 507     "and using the potential energy. This also allows for an entropy",
 508     "estimate using G = H - T S, where H is the enthalpy (H = U + p V)",
 509     "and S entropy."
 510   };
 511   static bool bSum=FALSE,bGibbs=FALSE,bAll=FALSE,bFluct=FALSE;
 512   static bool bDp=FALSE,bMutot=FALSE;
 513   static int  skip=0,nmol=1,ndf=3;
 514   static real reftemp=300.0,ezero=0;
 515   t_pargs pa[] = {
 516     { "-G",   FALSE, etBOOL,  {&bGibbs},
 517       "Do a free energy estimate" },
 518     { "-Gtemp", FALSE, etREAL,{&reftemp},
 519       "Reference temperature for free energy calculation" },
 520     { "-zero", FALSE, etREAL, {&ezero},
 521       "Subtract a zero-point energy" },
 522     { "-sum",  FALSE, etBOOL, {&bSum},
 523       "Sum the energy terms selected rather than display them all" },
 524     { "-dp",   FALSE, etBOOL, {&bDp},
 525       "Print energies in high precision" },
 526     { "-mutot",FALSE, etBOOL, {&bMutot},
 527       "Compute the total dipole moment from the components" },
 528     { "-skip", FALSE, etINT,  {&skip},
 529       "Skip number of frames between data points" },
 530     { "-aver", FALSE, etBOOL, {&bAll},
 531       "Print also the X1,t and sigma1,t, only if only 1 energy is requested" },
 532     { "-nmol", FALSE, etINT,  {&nmol},
 533       "Number of molecules in your sample: the energies are divided by this number" },
 534     { "-ndf",  FALSE, etINT,  {&ndf},
 535       "Number of degrees of freedom per molecule. Necessary for calculating the heat capacity" },
 536     { "-fluc", FALSE, etBOOL, {&bFluct},
 537       "Calculate autocorrelation of energy fluctuations rather than energy itself" }
 538   };
 539   static char *drleg[] = {
 540     "Running average",
 541     "Instantaneous"
 542   };
 543   static char *setnm[] = {
 544     "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
 545     "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
 546     "Volume",  "Pressure"
 547   };
 548
 549   FILE       *out,*fp_pairs=NULL;
 550   FILE       **drout;
 551   int        fp;
 552   int        ndr;
 553   int        timecheck=0;
 554   t_topology top;
 555   t_inputrec ir;
 556   t_energy   *oldee,**ee;
 557   t_drblock  dr;
 558   int        teller,teller_disre,nre,step[2],oldstep;
 559   real       t[2],oldt;
 560   int        cur=0;
 561 #define NEXT (1-cur)
 562   real       *bounds,*violaver=NULL;
 563   int        *index,*pair;
 564   int        nbounds=0,npairs;
 565   bool       bDisRe,bDRAll,bStarted,bCont,bEDR,bVisco;
 566   double     sum,sumaver,sumt;
 567   real       **eneset=NULL,*time=NULL;
 568   int        *set=NULL,i,j,k,nset,sss,nenergy;
 569   char       **enm=NULL,**leg=NULL,**pairleg;
 570   char       **nms;
 571   t_filenm   fnm[] = {
 572     { efENX, "-f",    NULL,      ffOPTRD },
 573     { efTPX, "-s",    NULL,      ffOPTRD },
 574     { efXVG, "-o",    "energy",  ffWRITE },
 575     { efXVG, "-viol", "violaver",ffOPTWR },
 576     { efXVG, "-pairs","pairs",   ffOPTWR },
 577     { efXVG, "-corr", "enecorr", ffOPTWR },
 578     { efXVG, "-vis",  "visco",   ffOPTWR }
 579   };
 580 #define NFILE asize(fnm)
 581   int     npargs;
 582   t_pargs *ppa;
 583
 584   CopyRight(stderr,argv[0]);
 585   npargs = asize(pa);
 586   ppa    = add_acf_pargs(&npargs,pa);
 587   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
 588                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL);
 589
 590   bDRAll = opt2bSet("-pairs",NFILE,fnm);
 591   bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
 592
 593   fp = open_enx(ftp2fn(efENX,NFILE,fnm),"r");
 594   do_enxnms(fp,&nre,&enm);
 595
 596   /* Initiate energies and set them to zero */
 597   dr.ndr=0;
 598   snew(ee,2);
 599   snew(ee[cur],nre);
 600   snew(ee[NEXT],nre);
 601   snew(oldee,nre);
 602   nenergy = 0;
 603
 604   bVisco = opt2bSet("-vis",NFILE,fnm);
 605
 606   if (!bDisRe) {
 607     if (bVisco) {
 608       nset=asize(setnm);
 609       snew(set,nset);
 610       /* This is nasty code... To extract Pres tensor, Volume and Temperature */
 611       for(j=0; (j<nset); j++) {
 612         for(i=0; (i<nre); i++) {
 613           if (strstr(enm[i],setnm[j])) {
 614             set[j]=i;
 615             break;
 616           }
 617         }
 618         if (i == nre)
 619           fatal_error(0,"Could not find term %s for viscosity calculation",
 620                       setnm[j]);
 621       }
 622     }
 623     else {
 624       set=select_it(nre,enm,&nset);
 625     }
 626     out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",
 627                  "E (kJ mol\\S-1\\N)");
 628
 629     snew(leg,nset+1);
 630     for(i=0; (i<nset); i++)
 631       leg[i]=enm[set[i]];
 632     if (bSum) {
 633       leg[nset]="Sum";
 634       xvgr_legend(out,nset+1,leg);
 635     }
 636     else
 637       xvgr_legend(out,nset,leg);
 638
 639     snew(eneset,nset+1);
 640     time = NULL;
 641   }
 642   else {
 643     nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
 644                        &top,&ir);
 645     snew(violaver,npairs);
 646     out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
 647                  "Time (ps)","nm");
 648     xvgr_legend(out,2,drleg);
 649     if (bDRAll) {
 650       fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
 651                         "Time (ps)","Distance (nm)");
 652       fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
 653               ir.dr_tau);
 654     }
 655   }
 656
 657   /* Initiate counters */
 658   teller       = 0;
 659   teller_disre = 0;
 660   step[cur]    = 0;
 661   t[cur]       = 0;
 662   oldstep      = -1;
 663   oldt         = 0;
 664   do {
 665     /* This loop searches for the first frame (when -b option is given),
 666      * or when this has been found it reads just one energy frame
 667      */
 668     do {
 669       bCont = do_enx(fp,&(t[NEXT]),&(step[NEXT]),&nre,ee[NEXT],&ndr,&dr);
 670
 671       if (bCont)
 672         timecheck = check_times(t[NEXT]);
 673
 674     } while (bCont && (timecheck < 0));
 675
 676     if ((timecheck == 0) && bCont) {
 677       if (nre > 0) {
 678         /*
 679          * Only copy the values we just read, if we are within the time bounds
 680          * It is necessary for statistics to start counting from 1
 681          */
 682         cur  = NEXT;
 683         step[cur]++;
 684
 685         if (oldstep == -1) {
 686           /*
 687            * Initiate the previous step data
 688            * Since step is incremented after reading, it is always >= 1,
 689            * therefore this will be executed only once.
 690            */
 691           oldstep = step[cur] - 1;
 692           oldt    = t[cur];
 693           /*
 694            * If we did not start at the first step, oldstep will be > 0
 695            * and we must initiate the data in the array. Otherwise it remains 0
 696            */
 697           if (oldstep > 0)
 698             for(i=0; (i<nre); i++) {
 699               oldee[i].esum = ee[cur][i].esum - ee[cur][i].e;
 700               oldee[i].eav  = ee[cur][i].eav  -
 701                 (sqr(oldee[i].esum - oldstep*ee[cur][i].e)/(oldstep*step[cur]));
 702             }
 703         }
 704       }
 705       /*
 706        * Define distance restraint legends. Can only be done after
 707        * the first frame has been read... (Then we know how many there are)
 708        */
 709       if (bDisRe && bDRAll && !leg && (ndr > 0)) {
 710         t_iatom   *fa;
 711         t_iparams *ip;
 712
 713         fa = top.idef.il[F_DISRES].iatoms;
 714         ip = top.idef.iparams;
 715
 716         if (ndr != top.idef.il[F_DISRES].nr/3)
 717           fatal_error(0,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
 718                       ndr,top.idef.il[F_DISRES].nr/3);
 719
 720         snew(pairleg,dr.ndr);
 721         for(i=0; (i<dr.ndr); i++) {
 722           snew(pairleg[i],30);
 723           j=fa[3*i+1];
 724           k=fa[3*i+2];
 725           sprintf(pairleg[i],"%d %s %d %s (%d)",
 726                   top.atoms.atom[j].resnr+1,*top.atoms.atomname[j],
 727                   top.atoms.atom[k].resnr+1,*top.atoms.atomname[k],
 728                   ip[fa[3*i]].disres.index);
 729         }
 730         set=select_it(dr.ndr,pairleg,&nset);
 731         snew(leg,2*nset);
 732         for(i=0; (i<nset); i++) {
 733           snew(leg[2*i],32);
 734           sprintf(leg[2*i],  "a %s",pairleg[set[i]]);
 735           snew(leg[2*i+1],32);
 736           sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
 737         }
 738         xvgr_legend(fp_pairs,2*nset,leg);
 739       }
 740
 741       /*
 742        * Store energies for analysis afterwards...
 743        */
 744       if (!bDisRe && (nre > 0)) {
 745         if ((nenergy % 1000) == 0) {
 746           srenew(time,nenergy+1000);
 747           for(i=0; (i<=nset); i++)
 748             srenew(eneset[i],nenergy+1000);
 749         }
 750         time[nenergy] = t[cur];
 751         sum=0;
 752           for(i=0; (i<nset); i++) {
 753             eneset[i][nenergy] = ee[cur][set[i]].e;
 754             sum += ee[cur][set[i]].e;
 755           }
 756         if (bSum)
 757           eneset[nset][nenergy] = sum;
 758         nenergy++;
 759       }
 760       /*
 761        * Printing time, only when we do not want to skip frames
 762        */
 763       if ((!skip) || ((teller % skip) == 0)) {
 764         if (bDisRe) {
 765           /*******************************************
 766            * D I S T A N C E   R E S T R A I N T S
 767            *******************************************/
 768           if (ndr > 0) {
 769             print_one(out,bDp,t[cur]);
 770             if (violaver == NULL)
 771               snew(violaver,dr.ndr);
 772
 773             /* Subtract bounds from distances, to calculate violations */
 774             calc_violations(dr.rt,dr.rav,
 775                             nbounds,pair,bounds,violaver,&sumt,&sumaver);
 776
 777             fprintf(out,"  %8.4f  %8.4f\n",sumaver,sumt);
 778             if (bDRAll) {
 779               print_one(fp_pairs,bDp,t[cur]);
 780               for(i=0; (i<nset); i++) {
 781                 sss=set[i];
 782                 fprintf(fp_pairs,"  %8.4f",mypow(dr.rav[sss],minthird));
 783                 fprintf(fp_pairs,"  %8.4f",dr.rt[sss]);
 784               }
 785               fprintf(fp_pairs,"\n");
 786             }
 787             teller_disre++;
 788           }
 789         }
 790         /*******************************************
 791          * E N E R G I E S
 792          *******************************************/
 793         else {
 794           if (nre > 0) {
 795             print_one(out,bDp,t[cur]);
 796             if (bSum)
 797               print_one(out,bDp,(eneset[nset][nenergy-1]-ezero)/nmol);
 798             else if ((nset == 1) && bAll) {
 799               print_one(out,bDp,ee[cur][set[0]].e);
 800               print_one(out,bDp,ee[cur][set[0]].esum);
 801               print_one(out,bDp,ee[cur][set[0]].eav);
 802             }
 803             else for(i=0; (i<nset); i++)
 804               print_one(out,bDp,(ee[cur][set[i]].e-ezero)/nmol);
 805
 806             fprintf(out,"\n");
 807           }
 808         }
 809       }
 810       teller++;
 811     }
 812   } while (bCont && (timecheck == 0));
 813
 814   fprintf(stderr,"\n");
 815   close_enx(fp);
 816
 817   ffclose(out);
 818   if (bDRAll)
 819     ffclose(fp_pairs);
 820
 821   if (bDisRe)
 822     analyse_disre(opt2fn("-viol",NFILE,fnm),
 823                   teller_disre,violaver,bounds,index,pair,nbounds);
 824   else
 825     analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
 826                  bGibbs,bSum,bFluct,bVisco,opt2fn("-vis",NFILE,fnm),
 827                  nmol,ndf,oldstep,oldt,step[cur],t[cur],time,reftemp,
 828                  oldee,ee[cur],nset,set,nenergy,eneset,leg);
 829
 830   xvgr_file(opt2fn("-o",NFILE,fnm),"-nxy");
 831
 832   thanx(stdout);
 833
 834   return 0;
 835 }