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 #include "mdebin.h"
  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 void chomp(char *buf)
 140 {
 141   int len = strlen(buf);
 142
 143   while ((len > 0) && (buf[len-1] == '\n')) {
 144     buf[len-1] = '\0';
 145     len--;
 146   }
 147 }
 148
 149 static int *select_by_name(int nre,gmx_enxnm_t *nm,int *nset)
 150 {
 151   gmx_bool *bE;
 152   int  n,k,kk,j,i,nmatch,nind,nss;
 153   int  *set;
 154   gmx_bool bEOF,bVerbose = TRUE,bLong=FALSE;
 155   char *ptr,buf[STRLEN];
 156   const char *fm4="%3d  %-14s";
 157   const char *fm2="%3d  %-34s";
 158   char **newnm=NULL;
 159
 160   if ((getenv("VERBOSE")) != NULL)
 161     bVerbose = FALSE;
 162
 163   fprintf(stderr,"\n");
 164   fprintf(stderr,"Select the terms you want from the following list by\n");
 165   fprintf(stderr,"selecting either (part of) the name or the number or a combination.\n");
 166   fprintf(stderr,"End your selection with an empty line or a zero.\n");
 167   fprintf(stderr,"-------------------------------------------------------------------\n");
 168
 169   snew(newnm,nre);
 170   j = 0;
 171   for(k=0; k<nre; k++) {
 172     newnm[k] = strdup(nm[k].name);
 173     /* Insert dashes in all the names */
 174     while ((ptr = strchr(newnm[k],' ')) != NULL) {
 175       *ptr = '-';
 176     }
 177     if ( bVerbose ) {
 178       if (j == 0) {
 179         if (k > 0) {
 180           fprintf(stderr,"\n");
 181         }
 182         bLong = FALSE;
 183         for(kk=k; kk<k+4; kk++) {
 184           if (kk < nre && strlen(nm[kk].name) > 14) {
 185             bLong = TRUE;
 186           }
 187         }
 188       } else {
 189         fprintf(stderr," ");
 190       }
 191       if (!bLong) {
 192         fprintf(stderr,fm4,k+1,newnm[k]);
 193         j++;
 194         if (j == 4) {
 195           j = 0;
 196         }
 197       } else {
 198         fprintf(stderr,fm2,k+1,newnm[k]);
 199         j++;
 200         if (j == 2) {
 201           j = 0;
 202         }
 203       }
 204     }
 205   }
 206   if ( bVerbose ) {
 207     fprintf(stderr,"\n\n");
 208   }
 209
 210   snew(bE,nre);
 211
 212   bEOF = FALSE;
 213   while (!bEOF && (fgets2(buf,STRLEN-1,stdin))) {
 214     /* Remove newlines */
 215     chomp(buf);
 216
 217     /* Remove spaces */
 218     trim(buf);
 219
 220     /* Empty line means end of input */
 221     bEOF = (strlen(buf) == 0);
 222     if (!bEOF) {
 223       ptr = buf;
 224       do {
 225         if (!bEOF) {
 226           /* First try to read an integer */
 227           nss   = sscanf(ptr,"%d",&nind);
 228           if (nss == 1) {
 229             /* Zero means end of input */
 230             if (nind == 0) {
 231               bEOF = TRUE;
 232             } else if ((1<=nind) && (nind<=nre)) {
 233               bE[nind-1] = TRUE;
 234             } else {
 235               fprintf(stderr,"number %d is out of range\n",nind);
 236             }
 237           }
 238           else {
 239             /* Now try to read a string */
 240             i = strlen(ptr);
 241             nmatch = 0;
 242             for(nind=0; nind<nre; nind++) {
 243               if (gmx_strcasecmp(newnm[nind],ptr) == 0) {
 244                 bE[nind] = TRUE;
 245                 nmatch++;
 246               }
 247             }
 248             if (nmatch == 0) {
 249               i = strlen(ptr);
 250               nmatch = 0;
 251               for(nind=0; nind<nre; nind++) {
 252                 if (gmx_strncasecmp(newnm[nind],ptr,i) == 0) {
 253                   bE[nind] = TRUE;
 254                   nmatch++;
 255                 }
 256               }
 257               if (nmatch == 0) {
 258                 fprintf(stderr,"String '%s' does not match anything\n",ptr);
 259               }
 260             }
 261           }
 262         }
 263         /* Look for the first space, and remove spaces from there */
 264         if ((ptr = strchr(ptr,' ')) != NULL) {
 265           trim(ptr);
 266         }
 267       } while (!bEOF && (ptr && (strlen(ptr) > 0)));
 268     }
 269   }
 270
 271   snew(set,nre);
 272   for(i=(*nset)=0; (i<nre); i++)
 273     if (bE[i])
 274       set[(*nset)++]=i;
 275
 276   sfree(bE);
 277
 278   if (*nset == 0)
 279     gmx_fatal(FARGS,"No energy terms selected");
 280
 281   for(i=0; (i<nre); i++)
 282     sfree(newnm[i]);
 283   sfree(newnm);
 284
 285   return set;
 286 }
 287
 288 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
 289 {
 290     gmx_mtop_t mtop;
 291     int        natoms;
 292     t_iatom    *iatom;
 293     matrix     box;
 294
 295     /* all we need is the ir to be able to write the label */
 296     read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,&mtop);
 297 }
 298
 299 static void get_orires_parms(const char *topnm,
 300                                int *nor,int *nex,int **label,real **obs)
 301 {
 302   gmx_mtop_t mtop;
 303   gmx_localtop_t *top;
 304   t_inputrec ir;
 305   t_iparams  *ip;
 306   int        natoms,i;
 307   t_iatom    *iatom;
 308   int        nb;
 309   matrix     box;
 310
 311   read_tpx(topnm,&ir,box,&natoms,NULL,NULL,NULL,&mtop);
 312   top = gmx_mtop_generate_local_top(&mtop,&ir);
 313
 314   ip       = top->idef.iparams;
 315   iatom    = top->idef.il[F_ORIRES].iatoms;
 316
 317   /* Count how many distance restraint there are... */
 318   nb = top->idef.il[F_ORIRES].nr;
 319   if (nb == 0)
 320     gmx_fatal(FARGS,"No orientation restraints in topology!\n");
 321
 322   *nor = nb/3;
 323   *nex = 0;
 324   snew(*label,*nor);
 325   snew(*obs,*nor);
 326   for(i=0; i<nb; i+=3) {
 327     (*label)[i/3] = ip[iatom[i]].orires.label;
 328     (*obs)[i/3]   = ip[iatom[i]].orires.obs;
 329     if (ip[iatom[i]].orires.ex >= *nex)
 330       *nex = ip[iatom[i]].orires.ex+1;
 331   }
 332   fprintf(stderr,"Found %d orientation restraints with %d experiments",
 333           *nor,*nex);
 334 }
 335
 336 static int get_bounds(const char *topnm,
 337                       real **bounds,int **index,int **dr_pair,int *npairs,
 338                       gmx_mtop_t *mtop,gmx_localtop_t **ltop,t_inputrec *ir)
 339 {
 340   gmx_localtop_t *top;
 341   t_functype *functype;
 342   t_iparams  *ip;
 343   int        natoms,i,j,k,type,ftype,natom;
 344   t_ilist    *disres;
 345   t_iatom    *iatom;
 346   real       *b;
 347   int        *ind,*pair;
 348   int        nb,label1;
 349   matrix     box;
 350
 351   read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,mtop);
 352   snew(*ltop,1);
 353   top = gmx_mtop_generate_local_top(mtop,ir);
 354   *ltop = top;
 355
 356   functype = top->idef.functype;
 357   ip       = top->idef.iparams;
 358
 359   /* Count how many distance restraint there are... */
 360   nb=top->idef.il[F_DISRES].nr;
 361   if (nb == 0)
 362     gmx_fatal(FARGS,"No distance restraints in topology!\n");
 363
 364   /* Allocate memory */
 365   snew(b,nb);
 366   snew(ind,nb);
 367   snew(pair,nb+1);
 368
 369   /* Fill the bound array */
 370   nb=0;
 371   for(i=0; (i<top->idef.ntypes); i++) {
 372     ftype = functype[i];
 373     if (ftype == F_DISRES) {
 374
 375       label1 = ip[i].disres.label;
 376       b[nb]   = ip[i].disres.up1;
 377       ind[nb] = label1;
 378       nb++;
 379     }
 380   }
 381   *bounds = b;
 382
 383   /* Fill the index array */
 384   label1  = -1;
 385   disres  = &(top->idef.il[F_DISRES]);
 386   iatom   = disres->iatoms;
 387   for(i=j=k=0; (i<disres->nr); ) {
 388     type  = iatom[i];
 389     ftype = top->idef.functype[type];
 390     natom = interaction_function[ftype].nratoms+1;
 391     if (label1 != top->idef.iparams[type].disres.label) {
 392       pair[j] = k;
 393       label1  = top->idef.iparams[type].disres.label;
 394       j ++;
 395     }
 396     k++;
 397     i += natom;
 398   }
 399   pair[j]  = k;
 400   *npairs = k;
 401   if (j != nb)
 402     gmx_incons("get_bounds for distance restraints");
 403
 404   *index   = ind;
 405   *dr_pair = pair;
 406
 407   return nb;
 408 }
 409
 410 static void calc_violations(real rt[],real rav3[],int nb,int index[],
 411                             real bounds[],real *viol,double *st,double *sa)
 412 {
 413   const   real sixth=1.0/6.0;
 414   int     i,j;
 415   double  rsum,rav,sumaver,sumt;
 416
 417   sumaver = 0;
 418   sumt    = 0;
 419   for(i=0; (i<nb); i++) {
 420     rsum = 0.0;
 421     rav  = 0.0;
 422     for(j=index[i]; (j<index[i+1]); j++) {
 423       if (viol)
 424         viol[j] += mypow(rt[j],-3.0);
 425       rav     += sqr(rav3[j]);
 426       rsum    += mypow(rt[j],-6);
 427     }
 428     rsum    = max(0.0,mypow(rsum,-sixth)-bounds[i]);
 429     rav     = max(0.0,mypow(rav, -sixth)-bounds[i]);
 430
 431     sumt    += rsum;
 432     sumaver += rav;
 433   }
 434   *st = sumt;
 435   *sa = sumaver;
 436 }
 437
 438 static void analyse_disre(const char *voutfn,    int nframes,
 439                           real violaver[], real bounds[], int index[],
 440                           int pair[],      int nbounds,
 441                           const output_env_t oenv)
 442 {
 443     FILE   *vout;
 444     double sum,sumt,sumaver;
 445     int    i,j;
 446
 447     /* Subtract bounds from distances, to calculate violations */
 448     calc_violations(violaver,violaver,
 449             nbounds,pair,bounds,NULL,&sumt,&sumaver);
 450
 451 #ifdef DEBUG
 452     fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
 453             sumaver);
 454     fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",
 455             sumt);
 456 #endif
 457     vout=xvgropen(voutfn,"r\\S-3\\N average violations","DR Index","nm",
 458             oenv);
 459     sum  = 0.0;
 460     sumt = 0.0;
 461     for(i=0; (i<nbounds); i++) {
 462         /* Do ensemble averaging */
 463         sumaver = 0;
 464         for(j=pair[i]; (j<pair[i+1]); j++)
 465             sumaver += sqr(violaver[j]/nframes);
 466         sumaver = max(0.0,mypow(sumaver,minsixth)-bounds[i]);
 467
 468         sumt   += sumaver;
 469         sum     = max(sum,sumaver);
 470         fprintf(vout,"%10d  %10.5e\n",index[i],sumaver);
 471     }
 472 #ifdef DEBUG
 473     for(j=0; (j<dr.ndr); j++)
 474         fprintf(vout,"%10d  %10.5e\n",j,mypow(violaver[j]/nframes,minthird));
 475 #endif
 476     ffclose(vout);
 477
 478     fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
 479             sumt);
 480     fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",sum);
 481
 482     do_view(oenv,voutfn,"-graphtype bar");
 483 }
 484
 485 static void einstein_visco(const char *fn,const char *fni,int nsets,
 486                            int nframes,real **sum,
 487                            real V,real T,int nsteps,double time[],
 488                            const output_env_t oenv)
 489 {
 490     FILE *fp0,*fp1;
 491     double av[4],avold[4];
 492     double fac,dt,di;
 493     int  i,j,m,nf4;
 494
 495     if (nframes < 1)
 496         return;
 497
 498     dt  = (time[1]-time[0]);
 499     nf4 = nframes/4+1;
 500
 501     for(i=0; i<=nsets; i++)
 502         avold[i] = 0;
 503     fp0=xvgropen(fni,"Shear viscosity integral",
 504             "Time (ps)","(kg m\\S-1\\N s\\S-1\\N ps)",oenv);
 505     fp1=xvgropen(fn,"Shear viscosity using Einstein relation",
 506             "Time (ps)","(kg m\\S-1\\N s\\S-1\\N)",oenv);
 507     for(i=1; i<nf4; i++) {
 508         fac = dt*nframes/nsteps;
 509         for(m=0; m<=nsets; m++)
 510             av[m] = 0;
 511         for(j=0; j<nframes-i; j++) {
 512             for(m=0; m<nsets; m++) {
 513                 di   = sqr(fac*(sum[m][j+i]-sum[m][j]));
 514
 515                 av[m]     += di;
 516                 av[nsets] += di/nsets;
 517             }
 518         }
 519         /* Convert to SI for the viscosity */
 520         fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
 521         fprintf(fp0,"%10g",time[i]-time[0]);
 522         for(m=0; (m<=nsets); m++) {
 523             av[m] = fac*av[m];
 524             fprintf(fp0,"  %10g",av[m]);
 525         }
 526         fprintf(fp0,"\n");
 527         fprintf(fp1,"%10g",0.5*(time[i]+time[i-1])-time[0]);
 528         for(m=0; (m<=nsets); m++) {
 529             fprintf(fp1,"  %10g",(av[m]-avold[m])/dt);
 530             avold[m] = av[m];
 531         }
 532         fprintf(fp1,"\n");
 533     }
 534     ffclose(fp0);
 535     ffclose(fp1);
 536 }
 537
 538 typedef struct {
 539     gmx_large_int_t np;
 540     double sum;
 541     double sav;
 542     double sav2;
 543 } ee_sum_t;
 544
 545 typedef struct {
 546     int b;
 547     ee_sum_t sum;
 548     gmx_large_int_t nst;
 549     gmx_large_int_t nst_min;
 550 } ener_ee_t;
 551
 552 static void clear_ee_sum(ee_sum_t *ees)
 553 {
 554     ees->sav  = 0;
 555     ees->sav2 = 0;
 556     ees->np   = 0;
 557     ees->sum  = 0;
 558 }
 559
 560 static void add_ee_sum(ee_sum_t *ees,double sum,int np)
 561 {
 562     ees->np  += np;
 563     ees->sum += sum;
 564 }
 565
 566 static void add_ee_av(ee_sum_t *ees)
 567 {
 568     double av;
 569
 570     av = ees->sum/ees->np;
 571     ees->sav  += av;
 572     ees->sav2 += av*av;
 573     ees->np    = 0;
 574     ees->sum   = 0;
 575 }
 576
 577 static double calc_ee2(int nb,ee_sum_t *ees)
 578 {
 579     return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
 580 }
 581
 582 static void set_ee_av(ener_ee_t *eee)
 583 {
 584     if (debug)
 585     {
 586         char buf[STEPSTRSIZE];
 587         fprintf(debug,"Storing average for err.est.: %s steps\n",
 588                 gmx_step_str(eee->nst,buf));
 589     }
 590     add_ee_av(&eee->sum);
 591     eee->b++;
 592     if (eee->b == 1 || eee->nst < eee->nst_min)
 593     {
 594         eee->nst_min = eee->nst;
 595     }
 596     eee->nst = 0;
 597 }
 598
 599 static void calc_averages(int nset,enerdata_t *edat,int nbmin,int nbmax)
 600 {
 601     int  nb,i,f,nee;
 602     double sum,sum2,sump,see2;
 603     gmx_large_int_t steps,np,p,bound_nb;
 604     enerdat_t *ed;
 605     exactsum_t *es;
 606     gmx_bool bAllZero;
 607     double x,sx,sy,sxx,sxy;
 608     ener_ee_t *eee;
 609
 610     /* Check if we have exact statistics over all points */
 611     for(i=0; i<nset; i++)
 612     {
 613         ed = &edat->s[i];
 614         ed->bExactStat = FALSE;
 615         if (edat->npoints > 0)
 616         {
 617             /* All energy file sum entries 0 signals no exact sums.
 618              * But if all energy values are 0, we still have exact sums.
 619              */
 620             bAllZero = TRUE;
 621             for(f=0; f<edat->nframes && !ed->bExactStat; f++)
 622             {
 623                 if (ed->ener[i] != 0)
 624                 {
 625                     bAllZero = FALSE;
 626                 }
 627                 ed->bExactStat = (ed->es[f].sum != 0);
 628             }
 629             if (bAllZero)
 630             {
 631                 ed->bExactStat = TRUE;
 632             }
 633         }
 634     }
 635
 636     snew(eee,nbmax+1);
 637     for(i=0; i<nset; i++)
 638     {
 639         ed = &edat->s[i];
 640
 641         sum  = 0;
 642         sum2 = 0;
 643         np   = 0;
 644         sx   = 0;
 645         sy   = 0;
 646         sxx  = 0;
 647         sxy  = 0;
 648         for(nb=nbmin; nb<=nbmax; nb++)
 649         {
 650             eee[nb].b     = 0;
 651             clear_ee_sum(&eee[nb].sum);
 652             eee[nb].nst   = 0;
 653             eee[nb].nst_min = 0;
 654         }
 655         for(f=0; f<edat->nframes; f++)
 656         {
 657             es = &ed->es[f];
 658
 659             if (ed->bExactStat)
 660             {
 661                 /* Add the sum and the sum of variances to the totals. */
 662                 p     = edat->points[f];
 663                 sump  = es->sum;
 664                 sum2 += es->sum2;
 665                 if (np > 0)
 666                 {
 667                     sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
 668                         *np*(np + p)/p;
 669                 }
 670             }
 671             else
 672             {
 673                 /* Add a single value to the sum and sum of squares. */
 674                 p     = 1;
 675                 sump  = ed->ener[f];
 676                 sum2 += dsqr(sump);
 677             }
 678
 679             /* sum has to be increased after sum2 */
 680             np  += p;
 681             sum += sump;
 682
 683             /* For the linear regression use variance 1/p.
 684              * Note that sump is the sum, not the average, so we don't need p*.
 685              */
 686             x    = edat->step[f] - 0.5*(edat->steps[f] - 1);
 687             sx  += p*x;
 688             sy  += sump;
 689             sxx += p*x*x;
 690             sxy += x*sump;
 691
 692             for(nb=nbmin; nb<=nbmax; nb++)
 693             {
 694                 /* Check if the current end step is closer to the desired
 695                  * block boundary than the next end step.
 696                  */
 697                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 698                 if (eee[nb].nst > 0 &&
 699                     bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
 700                 {
 701                     set_ee_av(&eee[nb]);
 702                 }
 703                 if (f == 0)
 704                 {
 705                     eee[nb].nst = 1;
 706                 }
 707                 else
 708                 {
 709                     eee[nb].nst += edat->step[f] - edat->step[f-1];
 710                 }
 711                 if (ed->bExactStat)
 712                 {
 713                     add_ee_sum(&eee[nb].sum,es->sum,edat->points[f]);
 714                 }
 715                 else
 716                 {
 717                     add_ee_sum(&eee[nb].sum,edat->s[i].ener[f],1);
 718                 }
 719                 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
 720                 if (edat->step[f]*nb >= bound_nb)
 721                 {
 722                     set_ee_av(&eee[nb]);
 723                 }
 724             }
 725         }
 726
 727         edat->s[i].av = sum/np;
 728         if (ed->bExactStat)
 729         {
 730             edat->s[i].rmsd = sqrt(sum2/np);
 731         }
 732         else
 733         {
 734             edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
 735         }
 736
 737         if (edat->nframes > 1)
 738         {
 739             edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
 740         }
 741         else
 742         {
 743             edat->s[i].slope = 0;
 744         }
 745
 746         nee  = 0;
 747         see2 = 0;
 748         for(nb=nbmin; nb<=nbmax; nb++)
 749         {
 750             /* Check if we actually got nb blocks and if the smallest
 751              * block is not shorter than 80% of the average.
 752              */
 753             if (debug)
 754             {
 755                 char buf1[STEPSTRSIZE],buf2[STEPSTRSIZE];
 756                 fprintf(debug,"Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
 757                         nb,eee[nb].b,
 758                         gmx_step_str(eee[nb].nst_min,buf1),
 759                         gmx_step_str(edat->nsteps,buf2));
 760             }
 761             if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
 762             {
 763                 see2 += calc_ee2(nb,&eee[nb].sum);
 764                 nee++;
 765             }
 766         }
 767         if (nee > 0)
 768         {
 769             edat->s[i].ee = sqrt(see2/nee);
 770         }
 771         else
 772         {
 773             edat->s[i].ee = -1;
 774         }
 775     }
 776     sfree(eee);
 777 }
 778
 779 static enerdata_t *calc_sum(int nset,enerdata_t *edat,int nbmin,int nbmax)
 780 {
 781     enerdata_t *esum;
 782     enerdat_t *s;
 783     int f,i;
 784     double sum;
 785
 786     snew(esum,1);
 787     *esum = *edat;
 788     snew(esum->s,1);
 789     s = &esum->s[0];
 790     snew(s->ener,esum->nframes);
 791     snew(s->es  ,esum->nframes);
 792
 793     s->bExactStat = TRUE;
 794     s->slope      = 0;
 795     for(i=0; i<nset; i++)
 796     {
 797         if (!edat->s[i].bExactStat)
 798         {
 799             s->bExactStat = FALSE;
 800         }
 801         s->slope += edat->s[i].slope;
 802     }
 803
 804     for(f=0; f<edat->nframes; f++)
 805     {
 806         sum = 0;
 807         for(i=0; i<nset; i++)
 808         {
 809             sum += edat->s[i].ener[f];
 810         }
 811         s->ener[f] = sum;
 812         sum = 0;
 813         for(i=0; i<nset; i++)
 814         {
 815             sum += edat->s[i].es[f].sum;
 816         }
 817         s->es[f].sum  = sum;
 818         s->es[f].sum2 = 0;
 819     }
 820
 821     calc_averages(1,esum,nbmin,nbmax);
 822
 823     return esum;
 824 }
 825
 826 static char *ee_pr(double ee,char *buf)
 827 {
 828     char   tmp[100];
 829     double rnd;
 830
 831     if (ee < 0)
 832     {
 833         sprintf(buf,"%s","--");
 834     }
 835     else
 836     {
 837         /* Round to two decimals by printing. */
 838         sprintf(tmp,"%.1e",ee);
 839         sscanf(tmp,"%lf",&rnd);
 840         sprintf(buf,"%g",rnd);
 841     }
 842
 843     return buf;
 844 }
 845
 846 static void remove_drift(int nset,int nbmin,int nbmax,real dt,enerdata_t *edat)
 847 {
 848 /* Remove the drift by performing a fit to y = ax+b.
 849    Uses 5 iterations. */
 850     int i,j,k;
 851     double delta,d,sd,sd2;
 852
 853     edat->npoints = edat->nframes;
 854     edat->nsteps = edat->nframes;
 855
 856     for(k=0; (k<5); k++)
 857     {
 858         for(i=0; (i<nset); i++)
 859         {
 860             delta = edat->s[i].slope*dt;
 861
 862             if (NULL != debug)
 863                 fprintf(debug,"slope for set %d is %g\n",i,edat->s[i].slope);
 864
 865             for(j=0; (j<edat->nframes); j++)
 866             {
 867                 edat->s[i].ener[j]   -= j*delta;
 868                 edat->s[i].es[j].sum  = 0;
 869                 edat->s[i].es[j].sum2 = 0;
 870             }
 871         }
 872         calc_averages(nset,edat,nbmin,nbmax);
 873     }
 874 }
 875
 876 static void calc_fluctuation_props(FILE *fp,
 877                                    gmx_bool bDriftCorr,real dt,
 878                                    int nset,int set[],int nmol,
 879                                    char **leg,enerdata_t *edat,
 880                                    int nbmin,int nbmax)
 881 {
 882     int i,j;
 883     double vvhh,vv,v,h,hh2,vv2,varv,hh,varh,tt,cv,cp,alpha,kappa,dcp,et,varet;
 884     double NANO3;
 885     enum { eVol, eEnth, eTemp, eEtot, eNR };
 886     const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
 887     int ii[eNR];
 888
 889     NANO3 = NANO*NANO*NANO;
 890     if (!bDriftCorr)
 891     {
 892         fprintf(fp,"\nYou may want to use the -driftcorr flag in order to correct\n"
 893                 "for spurious drift in the graphs. Note that this is not\n"
 894                 "a substitute for proper equilibration and sampling!\n");
 895     }
 896     else
 897     {
 898         remove_drift(nset,nbmin,nbmax,dt,edat);
 899     }
 900     for(i=0; (i<eNR); i++)
 901     {
 902         for(ii[i]=0; (ii[i]<nset &&
 903                       (gmx_strcasecmp(leg[ii[i]],my_ener[i]) != 0)); ii[i]++)
 904             ;
 905 /*        if (ii[i] < nset)
 906             fprintf(fp,"Found %s data.\n",my_ener[i]);
 907 */  }
 908     /* Compute it all! */
 909     vvhh = alpha = kappa = cp = dcp = cv = NOTSET;
 910
 911     /* Temperature */
 912     tt = NOTSET;
 913     if (ii[eTemp] < nset)
 914     {
 915         tt    = edat->s[ii[eTemp]].av;
 916     }
 917     /* Volume */
 918     vv = varv = NOTSET;
 919     if ((ii[eVol] < nset) && (ii[eTemp] < nset))
 920     {
 921         vv    = edat->s[ii[eVol]].av*NANO3;
 922         varv  = dsqr(edat->s[ii[eVol]].rmsd*NANO3);
 923         kappa = (varv/vv)/(BOLTZMANN*tt);
 924     }
 925     /* Enthalpy */
 926     hh = varh = NOTSET;
 927     if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
 928     {
 929         hh    = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
 930         varh  = dsqr(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
 931         cp    = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
 932     }
 933     /* Total energy */
 934     et = varet = NOTSET;
 935     if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
 936     {
 937         /* Only compute cv in constant volume runs, which we can test
 938            by checking whether the enthalpy was computed.
 939         */
 940         et    = edat->s[ii[eEtot]].av;
 941         varet = sqr(edat->s[ii[eEtot]].rmsd);
 942         cv    = KILO*((varet/nmol)/(BOLTZ*tt*tt));
 943     }
 944     /* Alpha, dcp */
 945     if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
 946     {
 947         vvhh = 0;
 948         for(j=0; (j<edat->nframes); j++)
 949         {
 950             v = edat->s[ii[eVol]].ener[j]*NANO3;
 951             h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
 952             vvhh += (v*h);
 953         }
 954         vvhh /= edat->nframes;
 955         alpha = (vvhh-vv*hh)/(vv*BOLTZMANN*tt*tt);
 956         dcp   = (vv*AVOGADRO/nmol)*tt*sqr(alpha)/(kappa);
 957     }
 958
 959     if (tt != NOTSET)
 960     {
 961         if (nmol < 2)
 962             fprintf(fp,"\nWARNING: nmol = %d, this may not be what you want.\n",
 963                     nmol);
 964         fprintf(fp,"\nTemperature dependent fluctuation properties at T = %g.\n",tt);
 965         fprintf(fp,"\nHeat capacities obtained from fluctuations do *not* include\n");
 966         fprintf(fp,"quantum corrections. If you want to get a more accurate estimate\n");
 967         fprintf(fp,"please use the g_dos program.\n\n");
 968         fprintf(fp,"WARNING: Please verify that your simulations are converged and perform\n"
 969                 "a block-averaging error analysis (not implemented in g_energy yet)\n");
 970
 971         if (debug != NULL)
 972         {
 973             if (varv != NOTSET)
 974                 fprintf(fp,"varv  =  %10g (m^6)\n",varv*AVOGADRO/nmol);
 975             if (vvhh != NOTSET)
 976                 fprintf(fp,"vvhh  =  %10g (m^3 J)\n",vvhh);
 977         }
 978         if (vv != NOTSET)
 979             fprintf(fp,"Volume                                   = %10g m^3/mol\n",
 980                     vv*AVOGADRO/nmol);
 981         if (varh != NOTSET)
 982             fprintf(fp,"Enthalpy                                 = %10g kJ/mol\n",
 983                 hh*AVOGADRO/(KILO*nmol));
 984         if (alpha != NOTSET)
 985             fprintf(fp,"Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
 986                     alpha);
 987         if (kappa != NOTSET)
 988         {
 989             fprintf(fp,"Isothermal Compressibility Kappa         = %10g (J/m^3)\n",
 990                     kappa);
 991             fprintf(fp,"Adiabatic bulk modulus                   = %10g (m^3/J)\n",
 992                     1.0/kappa);
 993         }
 994         if (cp != NOTSET)
 995             fprintf(fp,"Heat capacity at constant pressure Cp    = %10g J/mol K\n",
 996                     cp);
 997         if (cv != NOTSET)
 998             fprintf(fp,"Heat capacity at constant volume Cv      = %10g J/mol K\n",
 999                     cv);
1000         if (dcp != NOTSET)
1001             fprintf(fp,"Cp-Cv                                    =  %10g J/mol K\n",
1002                     dcp);
1003         please_cite(fp,"Allen1987a");
1004     }
1005     else
1006     {
1007         fprintf(fp,"You should select the temperature in order to obtain fluctuation properties.\n");
1008     }
1009 }
1010
1011 static void analyse_ener(gmx_bool bCorr,const char *corrfn,
1012                          gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,gmx_bool bTempFluct,
1013                          gmx_bool bVisco,const char *visfn,int nmol,
1014                          gmx_large_int_t start_step,double start_t,
1015                          gmx_large_int_t step,double t,
1016                          double time[], real reftemp,
1017                          enerdata_t *edat,
1018                          int nset,int set[],gmx_bool *bIsEner,
1019                          char **leg,gmx_enxnm_t *enm,
1020                          real Vaver,real ezero,
1021                          int nbmin,int nbmax,
1022                          const output_env_t oenv)
1023 {
1024   FILE *fp;
1025   /* Check out the printed manual for equations! */
1026   double Dt,aver,stddev,errest,delta_t,totaldrift;
1027   enerdata_t *esum=NULL;
1028   real xxx,integral,intBulk,Temp=0,Pres=0;
1029   real sfrac,oldfrac,diffsum,diffav,fstep,pr_aver,pr_stddev,pr_errest;
1030   double beta=0,expE,expEtot,*fee=NULL;
1031   gmx_large_int_t nsteps;
1032   int  nexact,nnotexact;
1033   double x1m,x1mk;
1034   int  i,j,k,nout;
1035   real chi2;
1036   char buf[256],eebuf[100];
1037
1038   nsteps  = step - start_step + 1;
1039   if (nsteps < 1) {
1040     fprintf(stdout,"Not enough steps (%s) for statistics\n",
1041             gmx_step_str(nsteps,buf));
1042   }
1043   else {
1044     /* Calculate the time difference */
1045     delta_t = t - start_t;
1046
1047     fprintf(stdout,"\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
1048             gmx_step_str(nsteps,buf),start_t,t,nset);
1049
1050     calc_averages(nset,edat,nbmin,nbmax);
1051
1052     if (bSum) {
1053         esum = calc_sum(nset,edat,nbmin,nbmax);
1054     }
1055
1056     if (edat->npoints == 0) {
1057       nexact    = 0;
1058       nnotexact = nset;
1059     } else {
1060       nexact    = 0;
1061       nnotexact = 0;
1062       for(i=0; (i<nset); i++) {
1063         if (edat->s[i].bExactStat) {
1064           nexact++;
1065         } else {
1066           nnotexact++;
1067         }
1068       }
1069     }
1070
1071     if (nnotexact == 0) {
1072       fprintf(stdout,"All statistics are over %s points\n",
1073               gmx_step_str(edat->npoints,buf));
1074     } else if (nexact == 0 || edat->npoints == edat->nframes) {
1075       fprintf(stdout,"All statistics are over %d points (frames)\n",
1076               edat->nframes);
1077     } else {
1078       fprintf(stdout,"The term%s",nnotexact==1 ? "" : "s");
1079       for(i=0; (i<nset); i++) {
1080         if (!edat->s[i].bExactStat) {
1081           fprintf(stdout," '%s'",leg[i]);
1082         }
1083       }
1084       fprintf(stdout," %s has statistics over %d points (frames)\n",
1085               nnotexact==1 ? "is" : "are",edat->nframes);
1086       fprintf(stdout,"All other statistics are over %s points\n",
1087               gmx_step_str(edat->npoints,buf));
1088     }
1089     fprintf(stdout,"\n");
1090
1091     fprintf(stdout,"%-24s %10s %10s %10s %10s",
1092             "Energy","Average","Err.Est.","RMSD","Tot-Drift");
1093     if (bFee)
1094       fprintf(stdout,"  %10s\n","-kT ln<e^(E/kT)>");
1095     else
1096       fprintf(stdout,"\n");
1097     fprintf(stdout,"-------------------------------------------------------------------------------\n");
1098
1099     /* Initiate locals, only used with -sum */
1100     expEtot=0;
1101     if (bFee) {
1102       beta = 1.0/(BOLTZ*reftemp);
1103       snew(fee,nset);
1104     }
1105     for(i=0; (i<nset); i++) {
1106       aver   = edat->s[i].av;
1107       stddev = edat->s[i].rmsd;
1108       errest = edat->s[i].ee;
1109
1110       if (bFee) {
1111         expE = 0;
1112         for(j=0; (j<edat->nframes); j++) {
1113           expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
1114         }
1115         if (bSum)
1116           expEtot+=expE/edat->nframes;
1117
1118         fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
1119       }
1120       if (strstr(leg[i],"empera") != NULL) {
1121         Temp = aver;
1122       } else if (strstr(leg[i],"olum") != NULL) {
1123         Vaver= aver;
1124       } else if (strstr(leg[i],"essure") != NULL) {
1125         Pres = aver;
1126       }
1127       if (bIsEner[i]) {
1128         pr_aver   = aver/nmol-ezero;
1129         pr_stddev = stddev/nmol;
1130         pr_errest = errest/nmol;
1131       }
1132       else {
1133         pr_aver   = aver;
1134         pr_stddev = stddev;
1135         pr_errest = errest;
1136       }
1137
1138       /* Multiply the slope in steps with the number of steps taken */
1139       totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1140       if (bIsEner[i])
1141       {
1142           totaldrift /= nmol;
1143       }
1144
1145       fprintf(stdout,"%-24s %10g %10s %10g %10g",
1146               leg[i],pr_aver,ee_pr(pr_errest,eebuf),pr_stddev,totaldrift);
1147       if (bFee)
1148         fprintf(stdout,"  %10g",fee[i]);
1149
1150       fprintf(stdout,"  (%s)\n",enm[set[i]].unit);
1151
1152       if (bFluct) {
1153         for(j=0; (j<edat->nframes); j++)
1154           edat->s[i].ener[j] -= aver;
1155       }
1156     }
1157     if (bSum) {
1158         totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1159       fprintf(stdout,"%-24s %10g %10s %10s %10g  (%s)",
1160               "Total",esum->s[0].av/nmol,ee_pr(esum->s[0].ee/nmol,eebuf),
1161               "--",totaldrift/nmol,enm[set[0]].unit);
1162       /* pr_aver,pr_stddev,a,totaldrift */
1163       if (bFee)
1164         fprintf(stdout,"  %10g  %10g\n",
1165                 log(expEtot)/beta + esum->s[0].av/nmol,log(expEtot)/beta);
1166       else
1167         fprintf(stdout,"\n");
1168     }
1169
1170     /* Do correlation function */
1171     if (edat->nframes > 1)
1172     {
1173         Dt = delta_t/(edat->nframes - 1);
1174     }
1175     else
1176     {
1177         Dt = 0;
1178     }
1179     if (bVisco) {
1180       const char* leg[] = { "Shear", "Bulk" };
1181       real factor;
1182       real **eneset;
1183       real **enesum;
1184
1185       /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1186
1187       /* Symmetrise tensor! (and store in first three elements)
1188        * And subtract average pressure!
1189        */
1190       snew(eneset,12);
1191       for(i=0; i<12; i++) {
1192           snew(eneset[i],edat->nframes);
1193       }
1194       snew(enesum,3);
1195       for(i=0; i<3; i++) {
1196         snew(enesum[i],edat->nframes);
1197       }
1198       for(i=0; (i<edat->nframes); i++) {
1199         eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1200         eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1201         eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1202         for(j=3; j<=11; j++) {
1203           eneset[j][i] = edat->s[j].ener[i];
1204         }
1205         eneset[11][i] -= Pres;
1206         enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1207         enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1208         enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1209       }
1210
1211       einstein_visco("evisco.xvg","eviscoi.xvg",
1212                      3,edat->nframes,enesum,Vaver,Temp,nsteps,time,oenv);
1213
1214       /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1215       /* Do it for shear viscosity */
1216       strcpy(buf,"Shear Viscosity");
1217       low_do_autocorr(corrfn,oenv,buf,edat->nframes,3,
1218                       (edat->nframes+1)/2,eneset,Dt,
1219                       eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1220
1221       /* Now for bulk viscosity */
1222       strcpy(buf,"Bulk Viscosity");
1223       low_do_autocorr(corrfn,oenv,buf,edat->nframes,1,
1224                       (edat->nframes+1)/2,&(eneset[11]),Dt,
1225                       eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1226
1227       factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1228       fp=xvgropen(visfn,buf,"Time (ps)","\\8h\\4 (cp)",oenv);
1229       xvgr_legend(fp,asize(leg),leg,oenv);
1230
1231       /* Use trapezium rule for integration */
1232       integral = 0;
1233       intBulk  = 0;
1234       nout = get_acfnout();
1235       if ((nout < 2) || (nout >= edat->nframes/2))
1236           nout = edat->nframes/2;
1237       for(i=1; (i<nout); i++)
1238       {
1239           integral += 0.5*(eneset[0][i-1]  + eneset[0][i])*factor;
1240           intBulk  += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1241           fprintf(fp,"%10g  %10g  %10g\n",(i*Dt),integral,intBulk);
1242       }
1243       ffclose(fp);
1244     }
1245     else if (bCorr) {
1246       if (bFluct)
1247           strcpy(buf,"Autocorrelation of Energy Fluctuations");
1248       else
1249           strcpy(buf,"Energy Autocorrelation");
1250 #if 0
1251       do_autocorr(corrfn,oenv,buf,edat->nframes,
1252                   bSum ? 1                 : nset,
1253                   bSum ? &edat->s[nset-1].ener : eneset,
1254                   (delta_t/edat->nframes),eacNormal,FALSE);
1255 #endif
1256     }
1257   }
1258 }
1259
1260 static void print_time(FILE *fp,double t)
1261 {
1262   fprintf(fp,"%12.6f",t);
1263 }
1264
1265 static void print1(FILE *fp,gmx_bool bDp,real e)
1266 {
1267   if (bDp)
1268     fprintf(fp,"  %16.12f",e);
1269   else
1270     fprintf(fp,"  %10.6f",e);
1271 }
1272
1273 static void fec(const char *ene2fn, const char *runavgfn,
1274                 real reftemp, int nset, int set[], char *leg[],
1275                 enerdata_t *edat, double time[],
1276                 const output_env_t oenv)
1277 {
1278   const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1279                            "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1280   FILE *fp;
1281   ener_file_t enx;
1282   int  nre,timecheck,step,nenergy,nenergy2,maxenergy;
1283   int  i,j;
1284   gmx_bool bCont;
1285   real aver, beta;
1286   real **eneset2;
1287   double dE, sum;
1288   gmx_enxnm_t *enm=NULL;
1289   t_enxframe *fr;
1290   char buf[22];
1291
1292   /* read second energy file */
1293   snew(fr,1);
1294   enm = NULL;
1295   enx = open_enx(ene2fn,"r");
1296   do_enxnms(enx,&(fr->nre),&enm);
1297
1298   snew(eneset2,nset+1);
1299   nenergy2=0;
1300   maxenergy=0;
1301   timecheck=0;
1302   do {
1303     /* This loop searches for the first frame (when -b option is given),
1304      * or when this has been found it reads just one energy frame
1305      */
1306     do {
1307       bCont = do_enx(enx,fr);
1308
1309       if (bCont)
1310         timecheck = check_times(fr->t);
1311
1312     } while (bCont && (timecheck < 0));
1313
1314     /* Store energies for analysis afterwards... */
1315     if ((timecheck == 0) && bCont) {
1316       if (fr->nre > 0) {
1317         if ( nenergy2 >= maxenergy ) {
1318           maxenergy += 1000;
1319           for(i=0; i<=nset; i++)
1320             srenew(eneset2[i],maxenergy);
1321         }
1322         if (fr->t != time[nenergy2])
1323           fprintf(stderr,"\nWARNING time mismatch %g!=%g at frame %s\n",
1324                   fr->t, time[nenergy2], gmx_step_str(fr->step,buf));
1325         for(i=0; i<nset; i++)
1326           eneset2[i][nenergy2] = fr->ener[set[i]].e;
1327         nenergy2++;
1328       }
1329     }
1330   } while (bCont && (timecheck == 0));
1331
1332   /* check */
1333   if (edat->nframes != nenergy2) {
1334     fprintf(stderr,"\nWARNING file length mismatch %d!=%d\n",
1335             edat->nframes,nenergy2);
1336   }
1337   nenergy = min(edat->nframes,nenergy2);
1338
1339   /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1340   fp=NULL;
1341   if (runavgfn) {
1342     fp=xvgropen(runavgfn,"Running average free energy difference",
1343                 "Time (" unit_time ")","\\8D\\4E (" unit_energy ")",oenv);
1344     xvgr_legend(fp,asize(ravgleg),ravgleg,oenv);
1345   }
1346   fprintf(stdout,"\n%-24s %10s\n",
1347           "Energy","dF = -kT ln < exp(-(EB-EA)/kT) >A");
1348   sum=0;
1349   beta = 1.0/(BOLTZ*reftemp);
1350   for(i=0; i<nset; i++) {
1351     if (gmx_strcasecmp(leg[i],enm[set[i]].name)!=0)
1352       fprintf(stderr,"\nWARNING energy set name mismatch %s!=%s\n",
1353               leg[i],enm[set[i]].name);
1354     for(j=0; j<nenergy; j++) {
1355       dE = eneset2[i][j] - edat->s[i].ener[j];
1356       sum += exp(-dE*beta);
1357       if (fp)
1358         fprintf(fp,"%10g %10g %10g\n",
1359                 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1360     }
1361     aver = -BOLTZ*reftemp*log(sum/nenergy);
1362     fprintf(stdout,"%-24s %10g\n",leg[i],aver);
1363   }
1364   if(fp) ffclose(fp);
1365   sfree(fr);
1366 }
1367
1368
1369 static void do_dhdl(t_enxframe *fr, t_inputrec *ir, FILE **fp_dhdl, const char *filename, gmx_bool bDp,
1370                     int *blocks, int *hists, int *samples, int *nlambdas, const output_env_t oenv)
1371 {
1372     const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
1373     char title[STRLEN],label_x[STRLEN],label_y[STRLEN], legend[STRLEN];
1374     char buf[STRLEN];
1375     gmx_bool first=FALSE;
1376     int nblock_hist=0, nblock_dh=0, nblock_dhcoll=0;
1377     int i,j,k;
1378     /* coll data */
1379     double temp=0, start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
1380     static int setnr=0;
1381     gmx_bool changing_lambda=FALSE;
1382
1383     /* now count the blocks & handle the global dh data */
1384     for(i=0;i<fr->nblock;i++)
1385     {
1386         if (fr->block[i].id == enxDHHIST)
1387         {
1388             nblock_hist++;
1389         }
1390         else if (fr->block[i].id == enxDH)
1391         {
1392             nblock_dh++;
1393         }
1394         else if (fr->block[i].id == enxDHCOLL)
1395         {
1396             nblock_dhcoll++;
1397             if ( (fr->block[i].nsub < 1) ||
1398                  (fr->block[i].sub[0].type != xdr_datatype_double) ||
1399                  (fr->block[i].sub[0].nr < 5))
1400             {
1401                 gmx_fatal(FARGS, "Unexpected block data");
1402             }
1403
1404             /* read the data from the DHCOLL block */
1405             temp =         fr->block[i].sub[0].dval[0];
1406             start_time =   fr->block[i].sub[0].dval[1];
1407             delta_time =   fr->block[i].sub[0].dval[2];
1408             start_lambda = fr->block[i].sub[0].dval[3];
1409             delta_lambda = fr->block[i].sub[0].dval[4];
1410             changing_lambda = (delta_lambda != 0);
1411         }
1412     }
1413
1414     if (nblock_hist == 0 && nblock_dh == 0)
1415     {
1416         /* don't do anything */
1417         return;
1418     }
1419     if (nblock_hist>0 && nblock_dh>0)
1420     {
1421         gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1422     }
1423     if (! *fp_dhdl )
1424     {
1425         if (nblock_dh>0)
1426         {
1427             /* we have standard, non-histogram data -- call open_dhdl to open the file */
1428             *fp_dhdl=open_dhdl(filename,ir,oenv);
1429         }
1430         else
1431         {
1432             sprintf(title,"N(%s)",deltag);
1433             sprintf(label_x,"%s (%s)",deltag,unit_energy);
1434             sprintf(label_y,"Samples");
1435             *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,oenv);
1436             sprintf(buf,"T = %g (K), %s = %g", temp, lambda, start_lambda);
1437             xvgr_subtitle(*fp_dhdl,buf,oenv);
1438         }
1439     }
1440
1441     (*hists)+=nblock_hist;
1442     (*blocks)+=nblock_dh;
1443     (*nlambdas) = nblock_hist+nblock_dh;
1444
1445     /* write the data */
1446     if (nblock_hist > 0)
1447     {
1448         gmx_large_int_t sum=0;
1449         /* histograms */
1450         for(i=0;i<fr->nblock;i++)
1451         {
1452             t_enxblock *blk=&(fr->block[i]);
1453             if (blk->id==enxDHHIST)
1454             {
1455                 double foreign_lambda, dx;
1456                 gmx_large_int_t x0;
1457                 int nhist, derivative;
1458
1459                 /* check the block types etc. */
1460                 if ( (blk->nsub < 2) ||
1461                      (blk->sub[0].type != xdr_datatype_double) ||
1462                      (blk->sub[1].type != xdr_datatype_large_int) ||
1463                      (blk->sub[0].nr < 2)  ||
1464                      (blk->sub[1].nr < 2) )
1465                 {
1466                     gmx_fatal(FARGS, "Unexpected block data in file");
1467                 }
1468                 foreign_lambda=blk->sub[0].dval[0];
1469                 dx=blk->sub[0].dval[1];
1470                 nhist=blk->sub[1].lval[0];
1471                 derivative=blk->sub[1].lval[1];
1472                 for(j=0;j<nhist;j++)
1473                 {
1474                     const char *lg[1];
1475                     x0=blk->sub[1].lval[2+j];
1476
1477                     if (!derivative)
1478                     {
1479                         sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1480                                 deltag, lambda, foreign_lambda,
1481                                 lambda, start_lambda);
1482                     }
1483                     else
1484                     {
1485                         sprintf(legend, "N(%s | %s=%g)",
1486                                 dhdl, lambda, start_lambda);
1487                     }
1488
1489                     lg[0]=legend;
1490                     xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1491                     setnr++;
1492                     for(k=0;k<blk->sub[j+2].nr;k++)
1493                     {
1494                         int hist;
1495                         double xmin, xmax;
1496
1497                         hist=blk->sub[j+2].ival[k];
1498                         xmin=(x0+k)*dx;
1499                         xmax=(x0+k+1)*dx;
1500                         fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
1501                                 xmax, hist);
1502                         sum+=hist;
1503                     }
1504                     /* multiple histogram data blocks in one histogram
1505                        mean that the second one is the reverse of the first one:
1506                        for dhdl derivatives, it's important to know both the
1507                        maximum and minimum values */
1508                     dx=-dx;
1509                 }
1510             }
1511         }
1512         (*samples) += (int)(sum/nblock_hist);
1513     }
1514     else
1515     {
1516         /* raw dh */
1517         int len=0;
1518         char **setnames=NULL;
1519         int nnames=nblock_dh;
1520
1521         for(i=0;i<fr->nblock;i++)
1522         {
1523             t_enxblock *blk=&(fr->block[i]);
1524             if (blk->id == enxDH)
1525             {
1526                 if (len == 0)
1527                 {
1528                     len=blk->sub[2].nr;
1529                 }
1530                 else
1531                 {
1532                     if (len!=blk->sub[2].nr)
1533                     {
1534                         gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1535                     }
1536                 }
1537             }
1538         }
1539         (*samples) += len;
1540
1541         for(i=0;i<len;i++)
1542         {
1543             double time=start_time + delta_time*i;
1544
1545             fprintf(*fp_dhdl,"%.4f ", time);
1546
1547             for(j=0;j<fr->nblock;j++)
1548             {
1549                 t_enxblock *blk=&(fr->block[j]);
1550                 if (blk->id == enxDH)
1551                 {
1552                     double value;
1553                     if (blk->sub[2].type == xdr_datatype_float)
1554                     {
1555                         value=blk->sub[2].fval[i];
1556                     }
1557                     else
1558                     {
1559                         value=blk->sub[2].dval[i];
1560                     }
1561                     /* we need to decide which data type it is based on the count*/
1562
1563                     if (j==1 && ir->bExpanded)
1564                     {
1565                         fprintf(*fp_dhdl,"%4d", (int)value);   /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
1566                     } else {
1567                         if (bDp) {
1568                             fprintf(*fp_dhdl," %#.12g", value);   /* print normal precision */
1569                         }
1570                         else
1571                         {
1572                             fprintf(*fp_dhdl," %#.8g", value);   /* print normal precision */
1573                         }
1574                     }
1575                 }
1576             }
1577             fprintf(*fp_dhdl, "\n");
1578         }
1579     }
1580 }
1581
1582
1583 int gmx_energy(int argc,char *argv[])
1584 {
1585   const char *desc[] = {
1586
1587     "[TT]g_energy[tt] extracts energy components or distance restraint",
1588     "data from an energy file. The user is prompted to interactively",
1589     "select the desired energy terms.[PAR]",
1590
1591     "Average, RMSD, and drift are calculated with full precision from the",
1592     "simulation (see printed manual). Drift is calculated by performing",
1593     "a least-squares fit of the data to a straight line. The reported total drift",
1594     "is the difference of the fit at the first and last point.",
1595     "An error estimate of the average is given based on a block averages",
1596     "over 5 blocks using the full-precision averages. The error estimate",
1597     "can be performed over multiple block lengths with the options",
1598     "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1599     "[BB]Note[bb] that in most cases the energy files contains averages over all",
1600     "MD steps, or over many more points than the number of frames in",
1601     "energy file. This makes the [TT]g_energy[tt] statistics output more accurate",
1602     "than the [TT].xvg[tt] output. When exact averages are not present in the energy",
1603     "file, the statistics mentioned above are simply over the single, per-frame",
1604     "energy values.[PAR]",
1605
1606     "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1607
1608     "Some fluctuation-dependent properties can be calculated provided",
1609     "the correct energy terms are selected, and that the command line option",
1610     "[TT]-fluct_props[tt] is given. The following properties",
1611     "will be computed:[BR]",
1612     "Property                        Energy terms needed[BR]",
1613     "---------------------------------------------------[BR]",
1614     "Heat capacity C[SUB]p[sub] (NPT sims):    Enthalpy, Temp     [BR]",
1615     "Heat capacity C[SUB]v[sub] (NVT sims):    Etot, Temp         [BR]",
1616     "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1617     "Isothermal compressibility:     Vol, Temp          [BR]",
1618     "Adiabatic bulk modulus:         Vol, Temp          [BR]",
1619     "---------------------------------------------------[BR]",
1620     "You always need to set the number of molecules [TT]-nmol[tt].",
1621     "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1622     "for quantum effects. Use the [TT]g_dos[tt] program if you need that (and you do).[PAR]"
1623     "When the [TT]-viol[tt] option is set, the time averaged",
1624     "violations are plotted and the running time-averaged and",
1625     "instantaneous sum of violations are recalculated. Additionally",
1626     "running time-averaged and instantaneous distances between",
1627     "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1628
1629     "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1630     "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1631     "The first two options plot the orientation, the last three the",
1632     "deviations of the orientations from the experimental values.",
1633     "The options that end on an 'a' plot the average over time",
1634     "as a function of restraint. The options that end on a 't'",
1635     "prompt the user for restraint label numbers and plot the data",
1636     "as a function of time. Option [TT]-odr[tt] plots the RMS",
1637     "deviation as a function of restraint.",
1638     "When the run used time or ensemble averaged orientation restraints,",
1639     "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1640     "not ensemble-averaged orientations and deviations instead of",
1641     "the time and ensemble averages.[PAR]",
1642
1643     "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1644     "tensor for each orientation restraint experiment. With option",
1645     "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1646
1647     "Option [TT]-odh[tt] extracts and plots the free energy data",
1648     "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1649     "from the [TT]ener.edr[tt] file.[PAR]",
1650
1651     "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1652     "difference with an ideal gas state: [BR]",
1653     "  [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
1654     "  [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln][BR]",
1655     "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1656     "the average is over the ensemble (or time in a trajectory).",
1657     "Note that this is in principle",
1658     "only correct when averaging over the whole (Boltzmann) ensemble",
1659     "and using the potential energy. This also allows for an entropy",
1660     "estimate using:[BR]",
1661     "  [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T[BR]",
1662     "  [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
1663     "[PAR]",
1664
1665     "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1666     "difference is calculated [BR] dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1667     "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1668     "files, and the average is over the ensemble A. The running average",
1669     "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1670     "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1671
1672   };
1673   static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE,bDriftCorr=FALSE;
1674   static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE,bFluctProps=FALSE;
1675   static int  skip=0,nmol=1,nbmin=5,nbmax=5;
1676   static real reftemp=300.0,ezero=0;
1677   t_pargs pa[] = {
1678     { "-fee",   FALSE, etBOOL,  {&bFee},
1679       "Do a free energy estimate" },
1680     { "-fetemp", FALSE, etREAL,{&reftemp},
1681       "Reference temperature for free energy calculation" },
1682     { "-zero", FALSE, etREAL, {&ezero},
1683       "Subtract a zero-point energy" },
1684     { "-sum",  FALSE, etBOOL, {&bSum},
1685       "Sum the energy terms selected rather than display them all" },
1686     { "-dp",   FALSE, etBOOL, {&bDp},
1687       "Print energies in high precision" },
1688     { "-nbmin", FALSE, etINT, {&nbmin},
1689       "Minimum number of blocks for error estimate" },
1690     { "-nbmax", FALSE, etINT, {&nbmax},
1691       "Maximum number of blocks for error estimate" },
1692     { "-mutot",FALSE, etBOOL, {&bMutot},
1693       "Compute the total dipole moment from the components" },
1694     { "-skip", FALSE, etINT,  {&skip},
1695       "Skip number of frames between data points" },
1696     { "-aver", FALSE, etBOOL, {&bPrAll},
1697       "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1698     { "-nmol", FALSE, etINT,  {&nmol},
1699       "Number of molecules in your sample: the energies are divided by this number" },
1700     { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1701       "Compute properties based on energy fluctuations, like heat capacity" },
1702     { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1703       "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1704     { "-fluc", FALSE, etBOOL, {&bFluct},
1705       "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1706     { "-orinst", FALSE, etBOOL, {&bOrinst},
1707       "Analyse instantaneous orientation data" },
1708     { "-ovec", FALSE, etBOOL, {&bOvec},
1709       "Also plot the eigenvectors with [TT]-oten[tt]" }
1710   };
1711   const char* drleg[] = {
1712     "Running average",
1713     "Instantaneous"
1714   };
1715   static const char *setnm[] = {
1716     "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1717     "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1718     "Volume",  "Pressure"
1719   };
1720
1721   FILE       *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1722   FILE       *fp_dhdl=NULL;
1723   FILE       **drout;
1724   ener_file_t fp;
1725   int        timecheck=0;
1726   gmx_mtop_t mtop;
1727   gmx_localtop_t *top=NULL;
1728   t_inputrec ir;
1729   t_energy   **ee;
1730   enerdata_t edat;
1731   gmx_enxnm_t *enm=NULL;
1732   t_enxframe *frame,*fr=NULL;
1733   int        cur=0;
1734 #define NEXT (1-cur)
1735   int        nre,teller,teller_disre,nfr;
1736   gmx_large_int_t start_step;
1737   int        nor=0,nex=0,norfr=0,enx_i=0;
1738   real       start_t;
1739   real       *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1740   int        *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1741   int        nbounds=0,npairs;
1742   gmx_bool       bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN,bDHDL;
1743   gmx_bool       bFoundStart,bCont,bEDR,bVisco;
1744   double     sum,sumaver,sumt,ener,dbl;
1745   double     *time=NULL;
1746   real       Vaver;
1747   int        *set=NULL,i,j,k,nset,sss;
1748   gmx_bool       *bIsEner=NULL;
1749   char       **pairleg,**odtleg,**otenleg;
1750   char       **leg=NULL;
1751   char       **nms;
1752   char       *anm_j,*anm_k,*resnm_j,*resnm_k;
1753   int        resnr_j,resnr_k;
1754   const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1755   char       buf[256];
1756   output_env_t oenv;
1757   t_enxblock *blk=NULL;
1758   t_enxblock *blk_disre=NULL;
1759   int        ndisre=0;
1760   int        dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas=0;
1761
1762   t_filenm   fnm[] = {
1763     { efEDR, "-f",    NULL,      ffREAD  },
1764     { efEDR, "-f2",   NULL,      ffOPTRD },
1765     { efTPX, "-s",    NULL,      ffOPTRD },
1766     { efXVG, "-o",    "energy",  ffWRITE },
1767     { efXVG, "-viol", "violaver",ffOPTWR },
1768     { efXVG, "-pairs","pairs",   ffOPTWR },
1769     { efXVG, "-ora",  "orienta", ffOPTWR },
1770     { efXVG, "-ort",  "orientt", ffOPTWR },
1771     { efXVG, "-oda",  "orideva", ffOPTWR },
1772     { efXVG, "-odr",  "oridevr", ffOPTWR },
1773     { efXVG, "-odt",  "oridevt", ffOPTWR },
1774     { efXVG, "-oten", "oriten",  ffOPTWR },
1775     { efXVG, "-corr", "enecorr", ffOPTWR },
1776     { efXVG, "-vis",  "visco",   ffOPTWR },
1777     { efXVG, "-ravg", "runavgdf",ffOPTWR },
1778     { efXVG, "-odh",  "dhdl"    ,ffOPTWR }
1779   };
1780 #define NFILE asize(fnm)
1781   int     npargs;
1782   t_pargs *ppa;
1783
1784   CopyRight(stderr,argv[0]);
1785   npargs = asize(pa);
1786   ppa    = add_acf_pargs(&npargs,pa);
1787   parse_common_args(&argc,argv,
1788                     PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1789                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1790
1791   bDRAll = opt2bSet("-pairs",NFILE,fnm);
1792   bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1793   bORA   = opt2bSet("-ora",NFILE,fnm);
1794   bORT   = opt2bSet("-ort",NFILE,fnm);
1795   bODA   = opt2bSet("-oda",NFILE,fnm);
1796   bODR   = opt2bSet("-odr",NFILE,fnm);
1797   bODT   = opt2bSet("-odt",NFILE,fnm);
1798   bORIRE = bORA || bORT || bODA || bODR || bODT;
1799   bOTEN  = opt2bSet("-oten",NFILE,fnm);
1800   bDHDL  = opt2bSet("-odh",NFILE,fnm);
1801
1802   nset = 0;
1803
1804   snew(frame,2);
1805   fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1806   do_enxnms(fp,&nre,&enm);
1807
1808   Vaver = -1;
1809
1810   bVisco = opt2bSet("-vis",NFILE,fnm);
1811
1812   if ((!bDisRe) && (!bDHDL))
1813   {
1814       if (bVisco) {
1815           nset=asize(setnm);
1816           snew(set,nset);
1817           /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1818           for(j=0; j<nset; j++) {
1819               for(i=0; i<nre; i++) {
1820                   if (strstr(enm[i].name,setnm[j])) {
1821                       set[j]=i;
1822                       break;
1823                   }
1824               }
1825               if (i == nre) {
1826                   if (gmx_strcasecmp(setnm[j],"Volume")==0) {
1827                       printf("Enter the box volume (" unit_volume "): ");
1828                       if(1 != scanf("%lf",&dbl))
1829                       {
1830                           gmx_fatal(FARGS,"Error reading user input");
1831                       }
1832                       Vaver = dbl;
1833                   } else
1834                       gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1835                                 setnm[j]);
1836               }
1837           }
1838       }
1839       else
1840       {
1841           set=select_by_name(nre,enm,&nset);
1842       }
1843       /* Print all the different units once */
1844       sprintf(buf,"(%s)",enm[set[0]].unit);
1845       for(i=1; i<nset; i++) {
1846           for(j=0; j<i; j++) {
1847               if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1848                   break;
1849               }
1850           }
1851           if (j == i) {
1852               strcat(buf,", (");
1853               strcat(buf,enm[set[i]].unit);
1854               strcat(buf,")");
1855           }
1856       }
1857       out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1858                    oenv);
1859
1860       snew(leg,nset+1);
1861       for(i=0; (i<nset); i++)
1862           leg[i] = enm[set[i]].name;
1863       if (bSum) {
1864           leg[nset]=strdup("Sum");
1865           xvgr_legend(out,nset+1,(const char**)leg,oenv);
1866       }
1867       else
1868           xvgr_legend(out,nset,(const char**)leg,oenv);
1869
1870       snew(bIsEner,nset);
1871       for(i=0; (i<nset); i++) {
1872           bIsEner[i] = FALSE;
1873           for (j=0; (j <= F_ETOT); j++)
1874               bIsEner[i] = bIsEner[i] ||
1875                         (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1876       }
1877
1878       if (bPrAll && nset > 1) {
1879           gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1880       }
1881
1882       time = NULL;
1883
1884       if (bORIRE || bOTEN)
1885           get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1886
1887       if (bORIRE) {
1888           if (bOrinst)
1889               enx_i = enxORI;
1890           else
1891               enx_i = enxOR;
1892
1893           if (bORA || bODA)
1894               snew(orient,nor);
1895           if (bODR)
1896               snew(odrms,nor);
1897           if (bORT || bODT) {
1898               fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1899               fprintf(stderr,"End your selection with 0\n");
1900               j = -1;
1901               orsel = NULL;
1902               do {
1903                   j++;
1904                   srenew(orsel,j+1);
1905                   if(1 != scanf("%d",&(orsel[j])))
1906                   {
1907                       gmx_fatal(FARGS,"Error reading user input");
1908                   }
1909               } while (orsel[j] > 0);
1910               if (orsel[0] == -1) {
1911                   fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1912                   norsel = nor;
1913                   srenew(orsel,nor);
1914                   for(i=0; i<nor; i++)
1915                       orsel[i] = i;
1916               } else {
1917                   /* Build the selection */
1918                   norsel=0;
1919                   for(i=0; i<j; i++) {
1920                       for(k=0; k<nor; k++)
1921                           if (or_label[k] == orsel[i]) {
1922                               orsel[norsel] = k;
1923                               norsel++;
1924                               break;
1925                           }
1926                       if (k == nor)
1927                           fprintf(stderr,"Orientation restraint label %d not found\n",
1928                                   orsel[i]);
1929                   }
1930               }
1931               snew(odtleg,norsel);
1932               for(i=0; i<norsel; i++) {
1933                   snew(odtleg[i],256);
1934                   sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1935               }
1936               if (bORT) {
1937                   fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1938                                 "Time (ps)","",oenv);
1939                   if (bOrinst)
1940                       fprintf(fort,"%s",orinst_sub);
1941                   xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
1942               }
1943               if (bODT) {
1944                   fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1945                                 "Orientation restraint deviation",
1946                                 "Time (ps)","",oenv);
1947                   if (bOrinst)
1948                       fprintf(fodt,"%s",orinst_sub);
1949                   xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
1950               }
1951           }
1952       }
1953       if (bOTEN) {
1954           foten=xvgropen(opt2fn("-oten",NFILE,fnm),
1955                          "Order tensor","Time (ps)","",oenv);
1956           snew(otenleg,bOvec ? nex*12 : nex*3);
1957           for(i=0; i<nex; i++) {
1958               for(j=0; j<3; j++) {
1959                   sprintf(buf,"eig%d",j+1);
1960                   otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
1961               }
1962               if (bOvec) {
1963                   for(j=0; j<9; j++) {
1964                       sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
1965                       otenleg[12*i+3+j] = strdup(buf);
1966                   }
1967               }
1968           }
1969           xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
1970       }
1971   }
1972   else if (bDisRe)
1973   {
1974       nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
1975                          &mtop,&top,&ir);
1976       snew(violaver,npairs);
1977       out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
1978                    "Time (ps)","nm",oenv);
1979       xvgr_legend(out,2,drleg,oenv);
1980       if (bDRAll) {
1981           fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
1982                             "Time (ps)","Distance (nm)",oenv);
1983           if (output_env_get_print_xvgr_codes(oenv))
1984               fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
1985                       ir.dr_tau);
1986       }
1987   } else if (bDHDL) {
1988       get_dhdl_parms(ftp2fn(efTPX,NFILE,fnm),&ir);
1989   }
1990
1991    /* Initiate energies and set them to zero */
1992    edat.nsteps  = 0;
1993    edat.npoints = 0;
1994    edat.nframes = 0;
1995    edat.step    = NULL;
1996    edat.steps   = NULL;
1997    edat.points  = NULL;
1998    snew(edat.s,nset);
1999
2000    /* Initiate counters */
2001    teller       = 0;
2002    teller_disre = 0;
2003    bFoundStart  = FALSE;
2004    start_step   = 0;
2005    start_t      = 0;
2006    do {
2007      /* This loop searches for the first frame (when -b option is given),
2008       * or when this has been found it reads just one energy frame
2009       */
2010      do {
2011          bCont = do_enx(fp,&(frame[NEXT]));
2012          if (bCont) {
2013              timecheck = check_times(frame[NEXT].t);
2014          }
2015      } while (bCont && (timecheck < 0));
2016
2017      if ((timecheck == 0) && bCont) {
2018        /* We read a valid frame, so we can use it */
2019        fr = &(frame[NEXT]);
2020
2021        if (fr->nre > 0) {
2022          /* The frame contains energies, so update cur */
2023          cur  = NEXT;
2024
2025              if (edat.nframes % 1000 == 0)
2026              {
2027                  srenew(edat.step,edat.nframes+1000);
2028                  memset(&(edat.step[edat.nframes]),0,1000*sizeof(edat.step[0]));
2029                  srenew(edat.steps,edat.nframes+1000);
2030                  memset(&(edat.steps[edat.nframes]),0,1000*sizeof(edat.steps[0]));
2031                  srenew(edat.points,edat.nframes+1000);
2032                  memset(&(edat.points[edat.nframes]),0,1000*sizeof(edat.points[0]));
2033
2034                  for(i=0; i<nset; i++)
2035                  {
2036                      srenew(edat.s[i].ener,edat.nframes+1000);
2037                      memset(&(edat.s[i].ener[edat.nframes]),0,
2038                             1000*sizeof(edat.s[i].ener[0]));
2039                      srenew(edat.s[i].es  ,edat.nframes+1000);
2040                      memset(&(edat.s[i].es[edat.nframes]),0,
2041                             1000*sizeof(edat.s[i].es[0]));
2042                  }
2043              }
2044
2045              nfr = edat.nframes;
2046              edat.step[nfr] = fr->step;
2047
2048              if (!bFoundStart)
2049              {
2050                  bFoundStart = TRUE;
2051                  /* Initiate the previous step data */
2052                  start_step = fr->step;
2053                  start_t    = fr->t;
2054                  /* Initiate the energy sums */
2055                  edat.steps[nfr]  = 1;
2056                  edat.points[nfr] = 1;
2057                  for(i=0; i<nset; i++)
2058                  {
2059                      sss = set[i];
2060                      edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2061                      edat.s[i].es[nfr].sum2 = 0;
2062                  }
2063                  edat.nsteps  = 1;
2064                  edat.npoints = 1;
2065              }
2066              else
2067              {
2068                  edat.steps[nfr] = fr->nsteps;
2069                  {
2070                      if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
2071                      {
2072                          if (fr->nsum <= 1)
2073                          {
2074                              edat.points[nfr] = 1;
2075                              for(i=0; i<nset; i++)
2076                              {
2077                                  sss = set[i];
2078                                  edat.s[i].es[nfr].sum  = fr->ener[sss].e;
2079                                  edat.s[i].es[nfr].sum2 = 0;
2080                              }
2081                              edat.npoints += 1;
2082                          }
2083                          else
2084                          {
2085                              edat.points[nfr] = fr->nsum;
2086                              for(i=0; i<nset; i++)
2087                              {
2088                                  sss = set[i];
2089                                  edat.s[i].es[nfr].sum  = fr->ener[sss].esum;
2090                                  edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2091                              }
2092                              edat.npoints += fr->nsum;
2093                          }
2094                      }
2095                      else
2096                      {
2097                          /* The interval does not match fr->nsteps:
2098                           * can not do exact averages.
2099                           */
2100                          edat.npoints = 0;
2101                      }
2102                      edat.nsteps = fr->step - start_step + 1;
2103                  }
2104              }
2105              for(i=0; i<nset; i++)
2106              {
2107                  edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2108              }
2109        }
2110        /*
2111         * Define distance restraint legends. Can only be done after
2112         * the first frame has been read... (Then we know how many there are)
2113         */
2114        blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
2115        if (bDisRe && bDRAll && !leg && blk_disre)
2116        {
2117            t_iatom   *fa;
2118            t_iparams *ip;
2119
2120            fa = top->idef.il[F_DISRES].iatoms;
2121            ip = top->idef.iparams;
2122            if (blk_disre->nsub != 2 ||
2123                (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2124            {
2125                gmx_incons("Number of disre sub-blocks not equal to 2");
2126            }
2127
2128            ndisre=blk_disre->sub[0].nr ;
2129            if (ndisre != top->idef.il[F_DISRES].nr/3)
2130            {
2131                gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2132                          ndisre,top->idef.il[F_DISRES].nr/3);
2133            }
2134            snew(pairleg,ndisre);
2135            for(i=0; i<ndisre; i++)
2136            {
2137                snew(pairleg[i],30);
2138                j=fa[3*i+1];
2139                k=fa[3*i+2];
2140                gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
2141                gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
2142                sprintf(pairleg[i],"%d %s %d %s (%d)",
2143                        resnr_j,anm_j,resnr_k,anm_k,
2144                        ip[fa[3*i]].disres.label);
2145            }
2146            set=select_it(ndisre,pairleg,&nset);
2147            snew(leg,2*nset);
2148            for(i=0; (i<nset); i++)
2149            {
2150                snew(leg[2*i],32);
2151                sprintf(leg[2*i],  "a %s",pairleg[set[i]]);
2152                snew(leg[2*i+1],32);
2153                sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
2154            }
2155            xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
2156        }
2157
2158        /*
2159         * Store energies for analysis afterwards...
2160         */
2161        if (!bDisRe && !bDHDL && (fr->nre > 0)) {
2162            if (edat.nframes % 1000 == 0) {
2163                srenew(time,edat.nframes+1000);
2164            }
2165            time[edat.nframes] = fr->t;
2166            edat.nframes++;
2167        }
2168        /*
2169         * Printing time, only when we do not want to skip frames
2170         */
2171        if (!skip || teller % skip == 0) {
2172      if (bDisRe) {
2173        /*******************************************
2174         * D I S T A N C E   R E S T R A I N T S
2175         *******************************************/
2176        if (ndisre > 0)
2177            {
2178  #ifndef GMX_DOUBLE
2179              float *disre_rt =     blk_disre->sub[0].fval;
2180              float *disre_rm3tav = blk_disre->sub[1].fval;
2181  #else
2182              double *disre_rt =     blk_disre->sub[0].dval;
2183              double *disre_rm3tav = blk_disre->sub[1].dval;
2184  #endif
2185
2186          print_time(out,fr->t);
2187          if (violaver == NULL)
2188            snew(violaver,ndisre);
2189
2190          /* Subtract bounds from distances, to calculate violations */
2191          calc_violations(disre_rt, disre_rm3tav,
2192                  nbounds,pair,bounds,violaver,&sumt,&sumaver);
2193
2194          fprintf(out,"  %8.4f  %8.4f\n",sumaver,sumt);
2195          if (bDRAll) {
2196            print_time(fp_pairs,fr->t);
2197            for(i=0; (i<nset); i++) {
2198          sss=set[i];
2199          fprintf(fp_pairs,"  %8.4f", mypow(disre_rm3tav[sss],minthird));
2200          fprintf(fp_pairs,"  %8.4f", disre_rt[sss]);
2201            }
2202            fprintf(fp_pairs,"\n");
2203          }
2204          teller_disre++;
2205        }
2206      }
2207      else if (bDHDL)
2208      {
2209          do_dhdl(fr, &ir, &fp_dhdl, opt2fn("-odh",NFILE,fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2210      }
2211
2212      /*******************************************
2213       * E N E R G I E S
2214       *******************************************/
2215      else {
2216          if (fr->nre > 0) {
2217              if (bPrAll)
2218              {
2219                  /* We skip frames with single points (usually only the first frame),
2220                   * since they would result in an average plot with outliers.
2221                   */
2222                  if (fr->nsum > 1) {
2223                      print_time(out,fr->t);
2224                       print1(out,bDp,fr->ener[set[0]].e);
2225                       print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
2226                       print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
2227                       fprintf(out,"\n");
2228                  }
2229              }
2230              else
2231              {
2232                  print_time(out,fr->t);
2233                  if (bSum)
2234                  {
2235                      sum = 0;
2236                      for(i=0; i<nset; i++)
2237                      {
2238                          sum += fr->ener[set[i]].e;
2239                      }
2240                      print1(out,bDp,sum/nmol-ezero);
2241                  }
2242                  else
2243                  {
2244                      for(i=0; (i<nset); i++)
2245                      {
2246                          if (bIsEner[i])
2247                          {
2248                              print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
2249                          }
2250                          else
2251                          {
2252                              print1(out,bDp,fr->ener[set[i]].e);
2253                          }
2254                      }
2255                  }
2256                  fprintf(out,"\n");
2257              }
2258          }
2259           blk = find_block_id_enxframe(fr, enx_i, NULL);
2260           if (bORIRE && blk)
2261           {
2262 #ifndef GMX_DOUBLE
2263               xdr_datatype dt=xdr_datatype_float;
2264 #else
2265               xdr_datatype dt=xdr_datatype_double;
2266 #endif
2267               real *vals;
2268
2269               if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2270               {
2271                   gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2272               }
2273 #ifndef GMX_DOUBLE
2274               vals=blk->sub[0].fval;
2275 #else
2276               vals=blk->sub[0].dval;
2277 #endif
2278
2279               if (blk->sub[0].nr != (size_t)nor)
2280                   gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2281               if (bORA || bODA)
2282               {
2283                   for(i=0; i<nor; i++)
2284                       orient[i] += vals[i];
2285               }
2286               if (bODR)
2287               {
2288                   for(i=0; i<nor; i++)
2289                       odrms[i] += sqr(vals[i]-oobs[i]);
2290               }
2291               if (bORT)
2292               {
2293                   fprintf(fort,"  %10f",fr->t);
2294                   for(i=0; i<norsel; i++)
2295                       fprintf(fort," %g",vals[orsel[i]]);
2296                   fprintf(fort,"\n");
2297               }
2298               if (bODT)
2299               {
2300                   fprintf(fodt,"  %10f",fr->t);
2301                   for(i=0; i<norsel; i++)
2302                       fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
2303                   fprintf(fodt,"\n");
2304               }
2305               norfr++;
2306           }
2307          blk = find_block_id_enxframe(fr, enxORT, NULL);
2308          if (bOTEN && blk)
2309          {
2310 #ifndef GMX_DOUBLE
2311              xdr_datatype dt=xdr_datatype_float;
2312 #else
2313              xdr_datatype dt=xdr_datatype_double;
2314 #endif
2315              real *vals;
2316
2317              if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2318                  gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2319 #ifndef GMX_DOUBLE
2320              vals=blk->sub[0].fval;
2321 #else
2322              vals=blk->sub[0].dval;
2323 #endif
2324
2325               if (blk->sub[0].nr != (size_t)(nex*12))
2326                   gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2327                             blk->sub[0].nr/12, nex);
2328               fprintf(foten,"  %10f",fr->t);
2329               for(i=0; i<nex; i++)
2330                   for(j=0; j<(bOvec?12:3); j++)
2331                       fprintf(foten," %g",vals[i*12+j]);
2332               fprintf(foten,"\n");
2333           }
2334         }
2335       }
2336       teller++;
2337     }
2338   } while (bCont && (timecheck == 0));
2339
2340   fprintf(stderr,"\n");
2341   close_enx(fp);
2342   if (out)
2343       ffclose(out);
2344
2345   if (bDRAll)
2346       ffclose(fp_pairs);
2347
2348   if (bORT)
2349       ffclose(fort);
2350   if (bODT)
2351       ffclose(fodt);
2352   if (bORA)
2353   {
2354       out = xvgropen(opt2fn("-ora",NFILE,fnm),
2355                      "Average calculated orientations",
2356                      "Restraint label","",oenv);
2357       if (bOrinst)
2358           fprintf(out,"%s",orinst_sub);
2359       for(i=0; i<nor; i++)
2360           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr);
2361       ffclose(out);
2362   }
2363   if (bODA) {
2364       out = xvgropen(opt2fn("-oda",NFILE,fnm),
2365                      "Average restraint deviation",
2366                      "Restraint label","",oenv);
2367       if (bOrinst)
2368           fprintf(out,"%s",orinst_sub);
2369       for(i=0; i<nor; i++)
2370           fprintf(out,"%5d  %g\n",or_label[i],orient[i]/norfr-oobs[i]);
2371       ffclose(out);
2372   }
2373   if (bODR) {
2374       out = xvgropen(opt2fn("-odr",NFILE,fnm),
2375                      "RMS orientation restraint deviations",
2376                      "Restraint label","",oenv);
2377       if (bOrinst)
2378           fprintf(out,"%s",orinst_sub);
2379       for(i=0; i<nor; i++)
2380           fprintf(out,"%5d  %g\n",or_label[i],sqrt(odrms[i]/norfr));
2381       ffclose(out);
2382   }
2383   if (bOTEN)
2384       ffclose(foten);
2385
2386   if (bDisRe)
2387   {
2388       analyse_disre(opt2fn("-viol",NFILE,fnm),
2389                     teller_disre,violaver,bounds,index,pair,nbounds,oenv);
2390   }
2391   else if (bDHDL)
2392   {
2393       if (fp_dhdl)
2394       {
2395           ffclose(fp_dhdl);
2396           printf("\n\nWrote %d lambda values with %d samples as ",
2397                  dh_lambdas, dh_samples);
2398           if (dh_hists > 0)
2399           {
2400               printf("%d dH histograms ", dh_hists);
2401           }
2402           if (dh_blocks> 0)
2403           {
2404               printf("%d dH data blocks ", dh_blocks);
2405           }
2406           printf("to %s\n", opt2fn("-odh",NFILE,fnm));
2407
2408       }
2409       else
2410       {
2411           gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f",NFILE,fnm));
2412       }
2413
2414   }
2415   else
2416   {
2417       double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2418       analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
2419                    bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
2420                    bVisco,opt2fn("-vis",NFILE,fnm),
2421                    nmol,
2422                    start_step,start_t,frame[cur].step,frame[cur].t,
2423                    time,reftemp,&edat,
2424                    nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
2425                    oenv);
2426       if (bFluctProps)
2427           calc_fluctuation_props(stdout,bDriftCorr,dt,nset,set,nmol,leg,&edat,
2428                                  nbmin,nbmax);
2429   }
2430   if (opt2bSet("-f2",NFILE,fnm)) {
2431       fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
2432           reftemp, nset, set, leg, &edat, time ,oenv);
2433   }
2434
2435   {
2436       const char *nxy = "-nxy";
2437
2438       do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
2439       do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
2440       do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
2441       do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
2442       do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
2443       do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
2444       do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
2445       do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
2446       do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
2447   }
2448   thanx(stderr);
2449
2450   return 0;
2451 }