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
1238 static void do_dhdl(t_enxframe *fr, FILE **fp_dhdl, const char *filename,
1239                     int *blocks, int *hists, int *samples, int *nlambdas,
1240                     const output_env_t oenv)
1241 {
1242     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
1243     char title[STRLEN],label_x[STRLEN],label_y[STRLEN], legend[STRLEN];
1244     char buf[STRLEN];
1245     gmx_bool first=FALSE;
1246     int nblock_hist=0, nblock_dh=0, nblock_dhcoll=0;
1247     int i,j,k;
1248     /* coll data */
1249     double temp=0, start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
1250     static int setnr=0;
1251
1252     /* now count the blocks & handle the global dh data */
1253     for(i=0;i<fr->nblock;i++)
1254     {
1255         if (fr->block[i].id == enxDHHIST)
1256         {
1257             nblock_hist++;
1258         }
1259         else if (fr->block[i].id == enxDH)
1260         {
1261             nblock_dh++;
1262         }
1263         else if (fr->block[i].id == enxDHCOLL)
1264         {
1265             nblock_dhcoll++;
1266             if ( (fr->block[i].nsub < 1) ||
1267                  (fr->block[i].sub[0].type != xdr_datatype_double) ||
1268                  (fr->block[i].sub[0].nr < 5))
1269             {
1270                 gmx_fatal(FARGS, "Unexpected block data");
1271             }
1272
1273             /* read the data from the DHCOLL block */
1274             temp =         fr->block[i].sub[0].dval[0];
1275             start_time =   fr->block[i].sub[0].dval[1];
1276             delta_time =   fr->block[i].sub[0].dval[2];
1277             start_lambda = fr->block[i].sub[0].dval[3];
1278             delta_lambda = fr->block[i].sub[0].dval[4];
1279         }
1280     }
1281
1282     if (nblock_hist == 0 && nblock_dh == 0)
1283     {
1284         /* don't do anything */
1285         return;
1286     }
1287     if (nblock_hist>0 && nblock_dh>0)
1288     {
1289         gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1290     }
1291     if (! *fp_dhdl )
1292     {
1293         if (nblock_dh>0)
1294         {
1295             sprintf(title,"%s, %s",dhdl,deltag);
1296             sprintf(label_x,"%s (%s)","Time",unit_time);
1297             sprintf(label_y,"(%s)",unit_energy);
1298         }
1299         else
1300         {
1301             sprintf(title,"N(%s)",deltag);
1302             sprintf(label_x,"%s (%s)",deltag,unit_energy);
1303             sprintf(label_y,"Samples");
1304         }
1305         *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,
1306                                oenv);
1307         sprintf(buf,"T = %g (K), %s = %g", temp, lambda, start_lambda);
1308         xvgr_subtitle(*fp_dhdl,buf,oenv);
1309         first=TRUE;
1310     }
1311
1312
1313
1314     (*hists)+=nblock_hist;
1315     (*blocks)+=nblock_dh;
1316     (*nlambdas) = nblock_hist+nblock_dh;
1317
1318
1319     /* write the data */
1320     if (nblock_hist > 0)
1321     {
1322         gmx_large_int_t sum=0;
1323         /* histograms */
1324         for(i=0;i<fr->nblock;i++)
1325         {
1326             t_enxblock *blk=&(fr->block[i]);
1327             if (blk->id==enxDHHIST)
1328             {
1329                 double foreign_lambda, dx;
1330                 gmx_large_int_t x0;
1331                 int nhist, derivative;
1332
1333                 /* check the block types etc. */
1334                 if ( (blk->nsub < 2) ||
1335                      (blk->sub[0].type != xdr_datatype_double) ||
1336                      (blk->sub[1].type != xdr_datatype_large_int) ||
1337                      (blk->sub[0].nr < 2)  ||
1338                      (blk->sub[1].nr < 2) )
1339                 {
1340                     gmx_fatal(FARGS, "Unexpected block data in file");
1341                 }
1342                 foreign_lambda=blk->sub[0].dval[0];
1343                 dx=blk->sub[0].dval[1];
1344                 nhist=blk->sub[1].lval[0];
1345                 derivative=blk->sub[1].lval[1];
1346                 for(j=0;j<nhist;j++)
1347                 {
1348                     const char *lg[1];
1349                     x0=blk->sub[1].lval[2+j];
1350
1351                     if (!derivative)
1352                     {
1353                         sprintf(legend,
1354                                 "N(%s(%s=%g) | %s=%g)",
1355                                 deltag, lambda, foreign_lambda,
1356                                 lambda, start_lambda);
1357                     }
1358                     else
1359                     {
1360                         sprintf(legend, "N(%s | %s=%g)",
1361                                 dhdl, lambda, start_lambda);
1362                     }
1363
1364                     lg[0]=legend;
1365                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1366                     setnr++;
1367                     for(k=0;k<blk->sub[j+2].nr;k++)
1368                     {
1369                         int hist;
1370                         double xmin, xmax;
1371
1372                         hist=blk->sub[j+2].ival[k];
1373                         xmin=(x0+k)*dx;
1374                         xmax=(x0+k+1)*dx;
1375                         fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
1376                                 xmax, hist);
1377                         sum+=hist;
1378                     }
1379                     /* multiple histogram data blocks in one histogram
1380                        mean that the second one is the reverse of the first one:
1381                        for dhdl derivatives, it's important to know both the
1382                        maximum and minimum values */
1383                     dx=-dx;
1384                 }
1385             }
1386         }
1387
1388         (*samples) += (int)(sum/nblock_hist);
1389     }
1390     else
1391     {
1392         /* raw dh */
1393         int len=0;
1394         char **setnames=NULL;
1395         int nnames=nblock_dh;
1396
1397         if (fabs(delta_lambda) > 1e-9)
1398         {
1399             nnames++;
1400         }
1401         if (first)
1402         {
1403             snew(setnames, nnames);
1404         }
1405         j=0;
1406
1407         if (fabs(delta_lambda) > 1e-9)
1408         {
1409             /* lambda is a plotted value */
1410             setnames[j]=gmx_strdup(lambda);
1411             j++;
1412         }
1413
1414
1415         for(i=0;i<fr->nblock;i++)
1416         {
1417             t_enxblock *blk=&(fr->block[i]);
1418             if (blk->id == enxDH)
1419             {
1420                 if (first)
1421                 {
1422                     /* do the legends */
1423                     int derivative;
1424                     double foreign_lambda;
1425
1426                     derivative=blk->sub[0].ival[0];
1427                     foreign_lambda=blk->sub[1].dval[0];
1428
1429                     if (derivative)
1430                     {
1431                         sprintf(buf, "%s %s %g",dhdl,lambda,start_lambda);
1432                     }
1433                     else
1434                     {
1435                         sprintf(buf, "%s %s %g",deltag,lambda,
1436                                 foreign_lambda);
1437                     }
1438                     setnames[j] = gmx_strdup(buf);
1439                     j++;
1440                 }
1441
1442                 if (len == 0)
1443                 {
1444                     len=blk->sub[2].nr;
1445                 }
1446                 else
1447                 {
1448                     if (len!=blk->sub[2].nr)
1449                     {
1450                         gmx_fatal(FARGS,
1451                                   "Length inconsistency in dhdl data");
1452                     }
1453                 }
1454             }
1455         }
1456
1457
1458         if (first)
1459         {
1460             xvgr_legend(*fp_dhdl, nblock_dh, (const char**)setnames, oenv);
1461             setnr += nblock_dh;
1462             for(i=0;i<nblock_dh;i++)
1463             {
1464                 sfree(setnames[i]);
1465             }
1466             sfree(setnames);
1467         }
1468
1469         (*samples) += len;
1470         for(i=0;i<len;i++)
1471         {
1472             double time=start_time + delta_time*i;
1473
1474             fprintf(*fp_dhdl,"%.4f", time);
1475             if (fabs(delta_lambda) > 1e-9)
1476             {
1477                 double lambda_now=i*delta_lambda + start_lambda;
1478                 fprintf(*fp_dhdl,"  %.4f", lambda_now);
1479             }
1480             for(j=0;j<fr->nblock;j++)
1481             {
1482                 t_enxblock *blk=&(fr->block[j]);
1483                 if (blk->id == enxDH)
1484                 {
1485                     double value;
1486                     if (blk->sub[2].type == xdr_datatype_float)
1487                     {
1488                         value=blk->sub[2].fval[i];
1489                     }
1490                     else
1491                     {
1492                         value=blk->sub[2].dval[i];
1493                     }
1494                     fprintf(*fp_dhdl,"  %g", value);
1495                 }
1496             }
1497             fprintf(*fp_dhdl, "\n");
1498         }
1499     }
1500 }
1501
1502
1503 int gmx_energy(int argc,char *argv[])
1504 {
1505   const char *desc[] = {
1506
1507     "[TT]g_energy[tt] extracts energy components or distance restraint",
1508     "data from an energy file. The user is prompted to interactively",
1509     "select the energy terms she wants.[PAR]",
1510
1511     "Average, RMSD and drift are calculated with full precision from the",
1512     "simulation (see printed manual). Drift is calculated by performing",
1513     "a least-squares fit of the data to a straight line. The reported total drift",
1514     "is the difference of the fit at the first and last point.",
1515     "An error estimate of the average is given based on a block averages",
1516     "over 5 blocks using the full precision averages. The error estimate",
1517     "can be performed over multiple block lengths with the options",
1518     "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1519     "Note that in most cases the energy files contains averages over all",
1520     "MD steps, or over many more points than the number of frames in",
1521     "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
1522     "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1523     "file the statistics mentioned above are simply over the single, per-frame",
1524     "energy values.[PAR]",
1525
1526     "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1527
1528     "Some fluctuation-dependent properties can be calculated provided",
1529     "the correct energy terms are selected. The following properties",
1530     "will be computed:[BR]",
1531     "Property                        Energy terms needed[BR]",
1532     "---------------------------------------------------[BR]",
1533     "Heat capacity Cp (NPT sims):    Enthalpy, Temp     [BR]",
1534     "Heat capacity Cv (NVT sims):    Etot, Temp         [BR]",
1535     "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1536     "Isothermal compressibility:     Vol, Temp          [BR]",
1537     "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1538     "---------------------------------------------------[BR]",
1539     "You always need to set the number of molecules [TT]-nmol[tt], and,",
1540     "if you used constraints in your simulations you will need to give",
1541     "the number of constraints per molecule [TT]-nconstr[tt] in order to",
1542     "correct for this: (nconstr/2) kB is subtracted from the heat",
1543     "capacity in this case. For instance in the case of rigid water",
1544     "you need to give the value 3 to this option.[PAR]",
1545
1546     "When the [TT]-viol[tt] option is set, the time averaged",
1547     "violations are plotted and the running time-averaged and",
1548     "instantaneous sum of violations are recalculated. Additionally",
1549     "running time-averaged and instantaneous distances between",
1550     "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1551
1552     "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1553     "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1554     "The first two options plot the orientation, the last three the",
1555     "deviations of the orientations from the experimental values.",
1556     "The options that end on an 'a' plot the average over time",
1557     "as a function of restraint. The options that end on a 't'",
1558     "prompt the user for restraint label numbers and plot the data",
1559     "as a function of time. Option [TT]-odr[tt] plots the RMS",
1560     "deviation as a function of restraint.",
1561     "When the run used time or ensemble averaged orientation restraints,",
1562     "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1563     "not ensemble-averaged orientations and deviations instead of",
1564     "the time and ensemble averages.[PAR]",
1565
1566     "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1567     "tensor for each orientation restraint experiment. With option",
1568     "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1569
1570     "Option [TT]-odh[tt] extracts and plots the free energy data",
1571     "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1572     "from the [TT]ener.edr[tt] file.[PAR]",
1573
1574     "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1575     "difference with an ideal gas state: [BR]",
1576     "  Delta A = A(N,V,T) - A_idgas(N,V,T) = kT ln < e^(Upot/kT) >[BR]",
1577     "  Delta G = G(N,p,T) - G_idgas(N,p,T) = kT ln < e^(Upot/kT) >[BR]",
1578     "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1579     "the average is over the ensemble (or time in a trajectory).",
1580     "Note that this is in principle",
1581     "only correct when averaging over the whole (Boltzmann) ensemble",
1582     "and using the potential energy. This also allows for an entropy",
1583     "estimate using:[BR]",
1584     "  Delta S(N,V,T) = S(N,V,T) - S_idgas(N,V,T) = (<Upot> - Delta A)/T[BR]",
1585     "  Delta S(N,p,T) = S(N,p,T) - S_idgas(N,p,T) = (<Upot> + pV - Delta G)/T",
1586     "[PAR]",
1587
1588     "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1589     "difference is calculated dF = -kT ln < e ^ -(EB-EA)/kT >A ,",
1590     "where EA and EB are the energies from the first and second energy",
1591     "files, and the average is over the ensemble A. [BB]NOTE[bb] that",
1592     "the energies must both be calculated from the same trajectory."
1593
1594   };
1595   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE;
1596   static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
1597   static int  skip=0,nmol=1,nconstr=0,nbmin=5,nbmax=5;
1598   static real reftemp=300.0,ezero=0;
1599   t_pargs pa[] = {
1600     { "-fee",   FALSE, etBOOL,  {&bFee},
1601       "Do a free energy estimate" },
1602     { "-fetemp", FALSE, etREAL,{&reftemp},
1603       "Reference temperature for free energy calculation" },
1604     { "-zero", FALSE, etREAL, {&ezero},
1605       "Subtract a zero-point energy" },
1606     { "-sum",  FALSE, etBOOL, {&bSum},
1607       "Sum the energy terms selected rather than display them all" },
1608     { "-dp",   FALSE, etBOOL, {&bDp},
1609       "Print energies in high precision" },
1610     { "-nbmin", FALSE, etINT, {&nbmin},
1611       "Minimum number of blocks for error estimate" },
1612     { "-nbmax", FALSE, etINT, {&nbmax},
1613       "Maximum number of blocks for error estimate" },
1614     { "-mutot",FALSE, etBOOL, {&bMutot},
1615       "Compute the total dipole moment from the components" },
1616     { "-skip", FALSE, etINT,  {&skip},
1617       "Skip number of frames between data points" },
1618     { "-aver", FALSE, etBOOL, {&bPrAll},
1619       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1620     { "-nmol", FALSE, etINT,  {&nmol},
1621       "Number of molecules in your sample: the energies are divided by this number" },
1622     { "-nconstr",  FALSE, etINT,  {&nconstr},
1623       "Number of constraints per molecule. Necessary for calculating the heat capacity" },
1624     { "-fluc", FALSE, etBOOL, {&bFluct},
1625       "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1626     { "-orinst", FALSE, etBOOL, {&bOrinst},
1627       "Analyse instantaneous orientation data" },
1628     { "-ovec", FALSE, etBOOL, {&bOvec},
1629       "Also plot the eigenvectors with [TT]-oten[tt]" }
1630   };
1631   const char* drleg[] = {
1632     "Running average",
1633     "Instantaneous"
1634   };
1635   static const char *setnm[] = {
1636     "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1637     "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1638     "Volume",  "Pressure"
1639   };
1640
1641   FILE       *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1642   FILE       *fp_dhdl=NULL;
1643   FILE       **drout;
1644   ener_file_t fp;
1645   int        timecheck=0;
1646   gmx_mtop_t mtop;
1647   gmx_localtop_t *top=NULL;
1648   t_inputrec ir;
1649   t_energy   **ee;
1650   enerdata_t edat;
1651   gmx_enxnm_t *enm=NULL;
1652   t_enxframe *frame,*fr=NULL;
1653   int        cur=0;
1654 #define NEXT (1-cur)
1655   int        nre,teller,teller_disre,nfr;
1656   gmx_large_int_t start_step;
1657   int        nor=0,nex=0,norfr=0,enx_i=0;
1658   real       start_t;
1659   real       *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1660   int        *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1661   int        nbounds=0,npairs;
1662   gmx_bool       bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN,bDHDL;
1663   gmx_bool       bFoundStart,bCont,bEDR,bVisco;
1664   double     sum,sumaver,sumt,ener,dbl;
1665   double     *time=NULL;
1666   real       Vaver;
1667   int        *set=NULL,i,j,k,nset,sss;
1668   gmx_bool       *bIsEner=NULL;
1669   char       **pairleg,**odtleg,**otenleg;
1670   char       **leg=NULL;
1671   char       **nms;
1672   char       *anm_j,*anm_k,*resnm_j,*resnm_k;
1673   int        resnr_j,resnr_k;
1674   const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1675   char       buf[256];
1676   output_env_t oenv;
1677   t_enxblock *blk=NULL;
1678   t_enxblock *blk_disre=NULL;
1679   int        ndisre=0;
1680   int        dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas=0;
1681
1682   t_filenm   fnm[] = {
1683     { efEDR, "-f",    NULL,      ffREAD  },
1684     { efEDR, "-f2",   NULL,      ffOPTRD },
1685     { efTPX, "-s",    NULL,      ffOPTRD },
1686     { efXVG, "-o",    "energy",  ffWRITE },
1687     { efXVG, "-viol", "violaver",ffOPTWR },
1688     { efXVG, "-pairs","pairs",   ffOPTWR },
1689     { efXVG, "-ora",  "orienta", ffOPTWR },
1690     { efXVG, "-ort",  "orientt", ffOPTWR },
1691     { efXVG, "-oda",  "orideva", ffOPTWR },
1692     { efXVG, "-odr",  "oridevr", ffOPTWR },
1693     { efXVG, "-odt",  "oridevt", ffOPTWR },
1694     { efXVG, "-oten", "oriten",  ffOPTWR },
1695     { efXVG, "-corr", "enecorr", ffOPTWR },
1696     { efXVG, "-vis",  "visco",   ffOPTWR },
1697     { efXVG, "-ravg", "runavgdf",ffOPTWR },
1698     { efXVG, "-odh",  "dhdl"    ,ffOPTWR }
1699   };
1700 #define NFILE asize(fnm)
1701   int     npargs;
1702   t_pargs *ppa;
1703
1704   CopyRight(stderr,argv[0]);
1705   npargs = asize(pa);
1706   ppa    = add_acf_pargs(&npargs,pa);
1707   parse_common_args(&argc,argv,
1708                     PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1709                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1710
1711   bDRAll = opt2bSet("-pairs",NFILE,fnm);
1712   bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1713   bORA   = opt2bSet("-ora",NFILE,fnm);
1714   bORT   = opt2bSet("-ort",NFILE,fnm);
1715   bODA   = opt2bSet("-oda",NFILE,fnm);
1716   bODR   = opt2bSet("-odr",NFILE,fnm);
1717   bODT   = opt2bSet("-odt",NFILE,fnm);
1718   bORIRE = bORA || bORT || bODA || bODR || bODT;
1719   bOTEN  = opt2bSet("-oten",NFILE,fnm);
1720   bDHDL  = opt2bSet("-odh",NFILE,fnm);
1721
1722   nset = 0;
1723
1724   snew(frame,2);
1725   fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1726   do_enxnms(fp,&nre,&enm);
1727
1728   Vaver = -1;
1729
1730   bVisco = opt2bSet("-vis",NFILE,fnm);
1731
1732   if (!bDisRe && !bDHDL)
1733   {
1734       if (bVisco) {
1735           nset=asize(setnm);
1736           snew(set,nset);
1737           /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1738           for(j=0; j<nset; j++) {
1739               for(i=0; i<nre; i++) {
1740                   if (strstr(enm[i].name,setnm[j])) {
1741                       set[j]=i;
1742                       break;
1743                   }
1744               }
1745               if (i == nre) {
1746                   if (gmx_strcasecmp(setnm[j],"Volume")==0) {
1747                       printf("Enter the box volume (" unit_volume "): ");
1748                       if(1 != scanf("%lf",&dbl))
1749                       {
1750                           gmx_fatal(FARGS,"Error reading user input");
1751                       }
1752                       Vaver = dbl;
1753                   } else
1754                       gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1755                                 setnm[j]);
1756               }
1757           }
1758       }
1759       else
1760       {
1761           set=select_by_name(nre,enm,&nset);
1762       }
1763       /* Print all the different units once */
1764       sprintf(buf,"(%s)",enm[set[0]].unit);
1765       for(i=1; i<nset; i++) {
1766           for(j=0; j<i; j++) {
1767               if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1768                   break;
1769               }
1770           }
1771           if (j == i) {
1772               strcat(buf,", (");
1773               strcat(buf,enm[set[i]].unit);
1774               strcat(buf,")");
1775           }
1776       }
1777       out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1778                    oenv);
1779
1780       snew(leg,nset+1);
1781       for(i=0; (i<nset); i++)
1782           leg[i] = enm[set[i]].name;
1783       if (bSum) {
1784           leg[nset]=strdup("Sum");
1785           xvgr_legend(out,nset+1,(const char**)leg,oenv);
1786       }
1787       else
1788           xvgr_legend(out,nset,(const char**)leg,oenv);
1789
1790       snew(bIsEner,nset);
1791       for(i=0; (i<nset); i++) {
1792           bIsEner[i] = FALSE;
1793           for (j=0; (j <= F_ETOT); j++)
1794               bIsEner[i] = bIsEner[i] ||
1795                         (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1796       }
1797
1798       if (bPrAll && nset > 1) {
1799           gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1800       }
1801
1802       time = NULL;
1803
1804       if (bORIRE || bOTEN)
1805           get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1806
1807       if (bORIRE) {
1808           if (bOrinst)
1809               enx_i = enxORI;
1810           else
1811               enx_i = enxOR;
1812
1813           if (bORA || bODA)
1814               snew(orient,nor);
1815           if (bODR)
1816               snew(odrms,nor);
1817           if (bORT || bODT) {
1818               fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1819               fprintf(stderr,"End your selection with 0\n");
1820               j = -1;
1821               orsel = NULL;
1822               do {
1823                   j++;
1824                   srenew(orsel,j+1);
1825                   if(1 != scanf("%d",&(orsel[j])))
1826                   {
1827                       gmx_fatal(FARGS,"Error reading user input");
1828                   }
1829               } while (orsel[j] > 0);
1830               if (orsel[0] == -1) {
1831                   fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1832                   norsel = nor;
1833                   srenew(orsel,nor);
1834                   for(i=0; i<nor; i++)
1835                       orsel[i] = i;
1836               } else {
1837                   /* Build the selection */
1838                   norsel=0;
1839                   for(i=0; i<j; i++) {
1840                       for(k=0; k<nor; k++)
1841                           if (or_label[k] == orsel[i]) {
1842                               orsel[norsel] = k;
1843                               norsel++;
1844                               break;
1845                           }
1846                       if (k == nor)
1847                           fprintf(stderr,"Orientation restraint label %d not found\n",
1848                                   orsel[i]);
1849                   }
1850               }
1851               snew(odtleg,norsel);
1852               for(i=0; i<norsel; i++) {
1853                   snew(odtleg[i],256);
1854                   sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1855               }
1856               if (bORT) {
1857                   fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1858                                 "Time (ps)","",oenv);
1859                   if (bOrinst)
1860                       fprintf(fort,"%s",orinst_sub);
1861                   xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
1862               }
1863               if (bODT) {
1864                   fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1865                                 "Orientation restraint deviation",
1866                                 "Time (ps)","",oenv);
1867                   if (bOrinst)
1868                       fprintf(fodt,"%s",orinst_sub);
1869                   xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
1870               }
1871           }
1872       }
1873       if (bOTEN) {
1874           foten=xvgropen(opt2fn("-oten",NFILE,fnm),
1875                          "Order tensor","Time (ps)","",oenv);
1876           snew(otenleg,bOvec ? nex*12 : nex*3);
1877           for(i=0; i<nex; i++) {
1878               for(j=0; j<3; j++) {
1879                   sprintf(buf,"eig%d",j+1);
1880                   otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
1881               }
1882               if (bOvec) {
1883                   for(j=0; j<9; j++) {
1884                       sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
1885                       otenleg[12*i+3+j] = strdup(buf);
1886                   }
1887               }
1888           }
1889           xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
1890       }
1891   }
1892   else if (bDisRe)
1893   {
1894       nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
1895                          &mtop,&top,&ir);
1896       snew(violaver,npairs);
1897       out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
1898                    "Time (ps)","nm",oenv);
1899       xvgr_legend(out,2,drleg,oenv);
1900       if (bDRAll) {
1901           fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
1902                             "Time (ps)","Distance (nm)",oenv);
1903           if (output_env_get_print_xvgr_codes(oenv))
1904               fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
1905                       ir.dr_tau);
1906       }
1907   }
1908
1909
1910   /* Initiate energies and set them to zero */
1911   edat.nsteps  = 0;
1912   edat.npoints = 0;
1913   edat.nframes = 0;
1914   edat.step    = NULL;
1915   edat.steps   = NULL;
1916   edat.points  = NULL;
1917   snew(edat.s,nset);
1918
1919   /* Initiate counters */
1920   teller       = 0;
1921   teller_disre = 0;
1922   bFoundStart  = FALSE;
1923   start_step   = 0;
1924   start_t      = 0;
1925   do {
1926     /* This loop searches for the first frame (when -b option is given),
1927      * or when this has been found it reads just one energy frame
1928      */
1929     do {
1930       bCont = do_enx(fp,&(frame[NEXT]));
1931
1932       if (bCont) {
1933         timecheck = check_times(frame[NEXT].t);
1934       }
1935     } while (bCont && (timecheck < 0));
1936
1937     if ((timecheck == 0) && bCont) {
1938       /* We read a valid frame, so we can use it */
1939       fr = &(frame[NEXT]);
1940
1941       if (fr->nre > 0) {
1942         /* The frame contains energies, so update cur */
1943         cur  = NEXT;
1944
1945                 if (edat.nframes % 1000 == 0)
1946             {
1947                 srenew(edat.step,edat.nframes+1000);
1948                 srenew(edat.steps,edat.nframes+1000);
1949                 srenew(edat.points,edat.nframes+1000);
1950                 for(i=0; i<nset; i++)
1951                 {
1952                     srenew(edat.s[i].ener,edat.nframes+1000);
1953                     srenew(edat.s[i].es  ,edat.nframes+1000);
1954                 }
1955             }
1956
1957                 nfr = edat.nframes;
1958             edat.step[nfr] = fr->step;
1959
1960             if (!bFoundStart)
1961             {
1962                 bFoundStart = TRUE;
1963                 /* Initiate the previous step data */
1964                 start_step = fr->step;
1965                 start_t    = fr->t;
1966                 /* Initiate the energy sums */
1967                 edat.steps[nfr]  = 1;
1968                 edat.points[nfr] = 1;
1969                 for(i=0; i<nset; i++)
1970                 {
1971                     sss = set[i];
1972                     edat.s[i].es[nfr].sum  = fr->ener[sss].e;
1973                     edat.s[i].es[nfr].sum2 = 0;
1974                 }
1975                 edat.nsteps  = 1;
1976                 edat.npoints = 1;
1977             }
1978             else
1979             {
1980                 edat.steps[nfr] = fr->nsteps;
1981                 {
1982                     if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
1983                     {
1984                         if (fr->nsum <= 1)
1985                         {
1986                             edat.points[nfr] = 1;
1987                             for(i=0; i<nset; i++)
1988                             {
1989                                 sss = set[i];
1990                                 edat.s[i].es[nfr].sum  = fr->ener[sss].e;
1991                                 edat.s[i].es[nfr].sum2 = 0;
1992                             }
1993                             edat.npoints += 1;
1994                         }
1995                         else
1996                         {
1997                             edat.points[nfr] = fr->nsum;
1998                             for(i=0; i<nset; i++)
1999                             {
2000                                 sss = set[i];
2001                                 edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2002                                 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2003                             }
2004                             edat.npoints += fr->nsum;
2005                         }
2006                     }
2007                     else
2008                     {
2009                         /* The interval does not match fr->nsteps:
2010                          * can not do exact averages.
2011                          */
2012                         edat.npoints = 0;
2013                     }
2014                     edat.nsteps = fr->step - start_step + 1;
2015                 }
2016             }
2017             for(i=0; i<nset; i++)
2018             {
2019                 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2020             }
2021       }
2022       /*
2023        * Define distance restraint legends. Can only be done after
2024        * the first frame has been read... (Then we know how many there are)
2025        */
2026       blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
2027       if (bDisRe && bDRAll && !leg && blk_disre)
2028       {
2029           t_iatom   *fa;
2030           t_iparams *ip;
2031
2032           fa = top->idef.il[F_DISRES].iatoms;
2033           ip = top->idef.iparams;
2034           if (blk_disre->nsub != 2 ||
2035               (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2036           {
2037               gmx_incons("Number of disre sub-blocks not equal to 2");
2038           }
2039
2040           ndisre=blk_disre->sub[0].nr ;
2041           if (ndisre != top->idef.il[F_DISRES].nr/3)
2042           {
2043               gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2044                         ndisre,top->idef.il[F_DISRES].nr/3);
2045           }
2046           snew(pairleg,ndisre);
2047           for(i=0; i<ndisre; i++)
2048           {
2049               snew(pairleg[i],30);
2050               j=fa[3*i+1];
2051               k=fa[3*i+2];
2052               gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
2053               gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
2054               sprintf(pairleg[i],"%d %s %d %s (%d)",
2055                       resnr_j,anm_j,resnr_k,anm_k,
2056                       ip[fa[3*i]].disres.label);
2057           }
2058           set=select_it(ndisre,pairleg,&nset);
2059           snew(leg,2*nset);
2060           for(i=0; (i<nset); i++)
2061           {
2062               snew(leg[2*i],32);
2063               sprintf(leg[2*i],  "a %s",pairleg[set[i]]);
2064               snew(leg[2*i+1],32);
2065               sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
2066           }
2067           xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
2068       }
2069
2070       /*
2071        * Store energies for analysis afterwards...
2072        */
2073       if (!bDisRe && !bDHDL && (fr->nre > 0)) {
2074         if (edat.nframes % 1000 == 0) {
2075           srenew(time,edat.nframes+1000);
2076         }
2077         time[edat.nframes] = fr->t;
2078         edat.nframes++;
2079       }
2080       /*
2081        * Printing time, only when we do not want to skip frames
2082        */
2083       if (!skip || teller % skip == 0) {
2084         if (bDisRe) {
2085           /*******************************************
2086            * D I S T A N C E   R E S T R A I N T S
2087            *******************************************/
2088           if (ndisre > 0)
2089           {
2090 #ifndef GMX_DOUBLE
2091             float *disre_rt =     blk_disre->sub[0].fval;
2092             float *disre_rm3tav = blk_disre->sub[1].fval;
2093 #else
2094             double *disre_rt =     blk_disre->sub[0].dval;
2095             double *disre_rm3tav = blk_disre->sub[1].dval;
2096 #endif
2097
2098             print_time(out,fr->t);
2099             if (violaver == NULL)
2100               snew(violaver,ndisre);
2101
2102             /* Subtract bounds from distances, to calculate violations */
2103             calc_violations(disre_rt, disre_rm3tav,
2104                             nbounds,pair,bounds,violaver,&sumt,&sumaver);
2105
2106             fprintf(out,"  %8.4f  %8.4f\n",sumaver,sumt);
2107             if (bDRAll) {
2108               print_time(fp_pairs,fr->t);
2109               for(i=0; (i<nset); i++) {
2110                 sss=set[i];
2111                 fprintf(fp_pairs,"  %8.4f", mypow(disre_rm3tav[sss],minthird));
2112                 fprintf(fp_pairs,"  %8.4f", disre_rt[sss]);
2113               }
2114               fprintf(fp_pairs,"\n");
2115             }
2116             teller_disre++;
2117           }
2118         }
2119         else if (bDHDL)
2120         {
2121             do_dhdl(fr, &fp_dhdl, opt2fn("-odh",NFILE,fnm),
2122                     &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas,
2123                     oenv);
2124         }
2125         /*******************************************
2126          * E N E R G I E S
2127          *******************************************/
2128         else {
2129           if (fr->nre > 0) {
2130             if (bPrAll)
2131             {
2132                 /* We skip frames with single points (usually only the first frame),
2133                  * since they would result in an average plot with outliers.
2134                  */
2135                 if (fr->nsum > 1) {
2136                     print_time(out,fr->t);
2137                      print1(out,bDp,fr->ener[set[0]].e);
2138                      print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
2139                      print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
2140                      fprintf(out,"\n");
2141                 }
2142             }
2143             else
2144             {
2145                 print_time(out,fr->t);
2146                 if (bSum)
2147                 {
2148                     sum = 0;
2149                     for(i=0; i<nset; i++)
2150                     {
2151                         sum += fr->ener[set[i]].e;
2152                     }
2153                     print1(out,bDp,sum/nmol-ezero);
2154                 }
2155                 else
2156                 {
2157                     for(i=0; (i<nset); i++)
2158                     {
2159                         if (bIsEner[i])
2160                         {
2161                             print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
2162                         }
2163                         else
2164                         {
2165                             print1(out,bDp,fr->ener[set[i]].e);
2166                         }
2167                     }
2168                 }
2169                 fprintf(out,"\n");
2170             }
2171           }
2172 #if 0
2173           /* we first count the blocks that have id 0: the orire blocks */
2174           block_orire=0;
2175           for(b=0;b<fr->nblock;b++)
2176           {
2177               if (fr->block[b].id == mde_block_type_orire)
2178                   nblock_orire++;
2179           }
2180 #endif
2181           blk = find_block_id_enxframe(fr, enx_i, NULL);
2182           if (bORIRE && blk)
2183           {
2184 #ifndef GMX_DOUBLE
2185               xdr_datatype dt=xdr_datatype_float;
2186 #else
2187               xdr_datatype dt=xdr_datatype_double;
2188 #endif
2189               real *vals;
2190
2191               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2192               {
2193                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2194               }
2195 #ifndef GMX_DOUBLE
2196               vals=blk->sub[0].fval;
2197 #else
2198               vals=blk->sub[0].dval;
2199 #endif
2200
2201               if (blk->sub[0].nr != (size_t)nor)
2202                   gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2203               if (bORA || bODA)
2204               {
2205                   for(i=0; i<nor; i++)
2206                       orient[i] += vals[i];
2207               }
2208               if (bODR)
2209               {
2210                   for(i=0; i<nor; i++)
2211                       odrms[i] += sqr(vals[i]-oobs[i]);
2212               }
2213               if (bORT)
2214               {
2215                   fprintf(fort,"  %10f",fr->t);
2216                   for(i=0; i<norsel; i++)
2217                       fprintf(fort," %g",vals[orsel[i]]);
2218                   fprintf(fort,"\n");
2219               }
2220               if (bODT)
2221               {
2222                   fprintf(fodt,"  %10f",fr->t);
2223                   for(i=0; i<norsel; i++)
2224                       fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
2225                   fprintf(fodt,"\n");
2226               }
2227               norfr++;
2228           }
2229           blk = find_block_id_enxframe(fr, enxORT, NULL);
2230           if (bOTEN && blk)
2231           {
2232 #ifndef GMX_DOUBLE
2233               xdr_datatype dt=xdr_datatype_float;
2234 #else
2235               xdr_datatype dt=xdr_datatype_double;
2236 #endif
2237               real *vals;
2238
2239               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2240                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2241 #ifndef GMX_DOUBLE
2242               vals=blk->sub[0].fval;
2243 #else
2244               vals=blk->sub[0].dval;
2245 #endif
2246
2247               if (blk->sub[0].nr != (size_t)(nex*12))
2248                   gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2249                             blk->sub[0].nr/12, nex);
2250               fprintf(foten,"  %10f",fr->t);
2251               for(i=0; i<nex; i++)
2252                   for(j=0; j<(bOvec?12:3); j++)
2253                       fprintf(foten," %g",vals[i*12+j]);
2254               fprintf(foten,"\n");
2255           }
2256         }
2257       }
2258     }
2259   } while (bCont && (timecheck == 0));
2260
2261   fprintf(stderr,"\n");
2262   close_enx(fp);
2263   if (out)
2264       ffclose(out);
2265
2266   if (bDRAll)
2267       ffclose(fp_pairs);
2268
2269   if (bORT)
2270       ffclose(fort);
2271   if (bODT)
2272       ffclose(fodt);
2273   if (bORA)
2274   {
2275       out = xvgropen(opt2fn("-ora",NFILE,fnm),
2276                      "Average calculated orientations",
2277                      "Restraint label","",oenv);
2278       if (bOrinst)
2279           fprintf(out,"%s",orinst_sub);
2280       for(i=0; i<nor; i++)
2281           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr);
2282       ffclose(out);
2283   }
2284   if (bODA) {
2285       out = xvgropen(opt2fn("-oda",NFILE,fnm),
2286                      "Average restraint deviation",
2287                      "Restraint label","",oenv);
2288       if (bOrinst)
2289           fprintf(out,"%s",orinst_sub);
2290       for(i=0; i<nor; i++)
2291           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr-oobs[i]);
2292       ffclose(out);
2293   }
2294   if (bODR) {
2295       out = xvgropen(opt2fn("-odr",NFILE,fnm),
2296                      "RMS orientation restraint deviations",
2297                      "Restraint label","",oenv);
2298       if (bOrinst)
2299           fprintf(out,"%s",orinst_sub);
2300       for(i=0; i<nor; i++)
2301           fprintf(out,"%5d  %g\n",or_label[i],sqrt(odrms[i]/norfr));
2302       ffclose(out);
2303   }
2304   if (bOTEN)
2305       ffclose(foten);
2306
2307   if (bDisRe)
2308   {
2309       analyse_disre(opt2fn("-viol",NFILE,fnm),
2310                     teller_disre,violaver,bounds,index,pair,nbounds,oenv);
2311   }
2312   else if (bDHDL)
2313   {
2314       if (fp_dhdl)
2315       {
2316           ffclose(fp_dhdl);
2317           printf("\n\nWrote %d lambda values with %d samples as ",
2318                  dh_lambdas, dh_samples);
2319           if (dh_hists > 0)
2320           {
2321               printf("%d dH histograms ", dh_hists);
2322           }
2323           if (dh_blocks> 0)
2324           {
2325               printf("%d dH data blocks ", dh_blocks);
2326           }
2327           printf("to %s\n", opt2fn("-odh",NFILE,fnm));
2328
2329       }
2330       else
2331       {
2332           gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f",NFILE,fnm));
2333       }
2334
2335   }
2336   else
2337   {
2338       analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
2339                    bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
2340                    bVisco,opt2fn("-vis",NFILE,fnm),
2341                    nmol,nconstr,start_step,start_t,frame[cur].step,frame[cur].t,
2342                    time,reftemp,&edat,
2343                    nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
2344                    oenv);
2345   }
2346   if (opt2bSet("-f2",NFILE,fnm)) {
2347       fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
2348           reftemp, nset, set, leg, &edat, time ,oenv);
2349   }
2350
2351   {
2352       const char *nxy = "-nxy";
2353
2354       do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
2355       do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
2356       do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
2357       do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
2358       do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
2359       do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
2360       do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
2361       do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
2362       do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
2363   }
2364   thanx(stderr);
2365
2366   return 0;
2367 }