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