src/tools/gmx_analyze.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 #include <math.h>
  40 #include <string.h>
  41 #include "statutil.h"
  42 #include "sysstuff.h"
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "macros.h"
  46 #include "gmx_fatal.h"
  47 #include "vec.h"
  48 #include "copyrite.h"
  49 #include "futil.h"
  50 #include "readinp.h"
  51 #include "statutil.h"
  52 #include "txtdump.h"
  53 #include "gstat.h"
  54 #include "gmx_matrix.h"
  55 #include "gmx_statistics.h"
  56 #include "xvgr.h"
  57 #include "gmx_ana.h"
  58 #include "geminate.h"
  59
  60 /* must correspond to char *avbar_opt[] declared in main() */
  61 enum { avbarSEL, avbarNONE, avbarSTDDEV, avbarERROR, avbar90, avbarNR };
  62
  63 static void power_fit(int n,int nset,real **val,real *t)
  64 {
  65   real *x,*y,quality,a,b,r;
  66   int  s,i;
  67
  68   snew(x,n);
  69   snew(y,n);
  70
  71   if (t[0]>0) {
  72     for(i=0; i<n; i++)
  73       if (t[0]>0)
  74         x[i] = log(t[i]);
  75   } else {
  76     fprintf(stdout,"First time is not larger than 0, using index number as time for power fit\n");
  77     for(i=0; i<n; i++)
  78       x[i] = log(i+1);
  79   }
  80
  81   for(s=0; s<nset; s++) {
  82     i=0;
  83     for(i=0; i<n && val[s][i]>=0; i++)
  84       y[i] = log(val[s][i]);
  85     if (i < n)
  86       fprintf(stdout,"Will power fit up to point %d, since it is not larger than 0\n",i);
  87     lsq_y_ax_b(i,x,y,&a,&b,&r,&quality);
  88     fprintf(stdout,"Power fit set %3d:  error %.3f  a %g  b %g\n",
  89             s+1,quality,a,exp(b));
  90   }
  91
  92   sfree(y);
  93   sfree(x);
  94 }
  95
  96 static real cosine_content(int nhp,int n,real *y)
  97      /* Assumes n equidistant points */
  98 {
  99   double fac,cosyint,yyint;
 100   int i;
 101
 102   if (n < 2)
 103     return 0;
 104
 105   fac = M_PI*nhp/(n-1);
 106
 107   cosyint = 0;
 108   yyint = 0;
 109   for(i=0; i<n; i++) {
 110     cosyint += cos(fac*i)*y[i];
 111     yyint += y[i]*y[i];
 112   }
 113
 114   return 2*cosyint*cosyint/(n*yyint);
 115 }
 116
 117 static void plot_coscont(const char *ccfile,int n,int nset,real **val,
 118                          const output_env_t oenv)
 119 {
 120   FILE *fp;
 121   int  s;
 122   real cc;
 123
 124   fp = xvgropen(ccfile,"Cosine content","set / half periods","cosine content",
 125                 oenv);
 126
 127   for(s=0; s<nset; s++) {
 128     cc = cosine_content(s+1,n,val[s]);
 129     fprintf(fp," %d %g\n",s+1,cc);
 130     fprintf(stdout,"Cosine content of set %d with %.1f periods: %g\n",
 131             s+1,0.5*(s+1),cc);
 132   }
 133   fprintf(stdout,"\n");
 134
 135   ffclose(fp);
 136 }
 137
 138 static void regression_analysis(int n,gmx_bool bXYdy,
 139                                 real *x,int nset,real **val)
 140 {
 141   real S,chi2,a,b,da,db,r=0;
 142   int ok;
 143
 144   if (bXYdy || (nset == 1))
 145   {
 146       printf("Fitting data to a function f(x) = ax + b\n");
 147       printf("Minimizing residual chi2 = Sum_i w_i [f(x_i) - y_i]2\n");
 148       printf("Error estimates will be given if w_i (sigma) values are given\n");
 149       printf("(use option -xydy).\n\n");
 150       if (bXYdy)
 151       {
 152           if ((ok = lsq_y_ax_b_error(n,x,val[0],val[1],&a,&b,&da,&db,&r,&S)) != estatsOK)
 153               gmx_fatal(FARGS,"Error fitting the data: %s",
 154                         gmx_stats_message(ok));
 155       }
 156       else
 157       {
 158           if ((ok = lsq_y_ax_b(n,x,val[0],&a,&b,&r,&S)) != estatsOK)
 159               gmx_fatal(FARGS,"Error fitting the data: %s",
 160                         gmx_stats_message(ok));
 161       }
 162       chi2 = sqr((n-2)*S);
 163       printf("Chi2                    = %g\n",chi2);
 164       printf("S (Sqrt(Chi2/(n-2))     = %g\n",S);
 165       printf("Correlation coefficient = %.1f%%\n",100*r);
 166       printf("\n");
 167       if (bXYdy) {
 168           printf("a    = %g +/- %g\n",a,da);
 169           printf("b    = %g +/- %g\n",b,db);
 170       }
 171       else {
 172           printf("a    = %g\n",a);
 173           printf("b    = %g\n",b);
 174       }
 175   }
 176   else
 177   {
 178       double chi2,*a,**xx,*y;
 179       int i,j;
 180
 181       snew(y,n);
 182       snew(xx,nset-1);
 183       for(j=0; (j<nset-1); j++)
 184           snew(xx[j],n);
 185       for(i=0; (i<n); i++)
 186       {
 187           y[i] = val[0][i];
 188           for(j=1; (j<nset); j++)
 189               xx[j-1][i] = val[j][i];
 190       }
 191       snew(a,nset-1);
 192       chi2 = multi_regression(NULL,n,y,nset-1,xx,a);
 193       printf("Fitting %d data points in %d sets\n",n,nset-1);
 194       printf("chi2 = %g\n",chi2);
 195       printf("A =");
 196       for(i=0; (i<nset-1); i++)
 197       {
 198           printf("  %g",a[i]);
 199           sfree(xx[i]);
 200       }
 201       printf("\n");
 202       sfree(xx);
 203       sfree(y);
 204       sfree(a);
 205   }
 206 }
 207
 208 void histogram(const char *distfile,real binwidth,int n, int nset, real **val,
 209                const output_env_t oenv)
 210 {
 211   FILE *fp;
 212   int  i,s;
 213   double min,max;
 214   int  nbin;
 215 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
 216   long long int *histo;
 217 #else
 218   double        *histo;
 219 #endif
 220
 221   min = val[0][0];
 222   max = val[0][0];
 223   for(s=0; s<nset; s++)
 224     for(i=0; i<n; i++)
 225       if (val[s][i] < min)
 226         min = val[s][i];
 227       else if (val[s][i] > max)
 228         max = val[s][i];
 229
 230   min = binwidth*floor(min/binwidth);
 231   max = binwidth*ceil(max/binwidth);
 232   if (min != 0)
 233     min -= binwidth;
 234   max += binwidth;
 235
 236   nbin = (int)((max - min)/binwidth + 0.5) + 1;
 237   fprintf(stderr,"Making distributions with %d bins\n",nbin);
 238   snew(histo,nbin);
 239   fp = xvgropen(distfile,"Distribution","","",oenv);
 240   for(s=0; s<nset; s++) {
 241     for(i=0; i<nbin; i++)
 242       histo[i] = 0;
 243     for(i=0; i<n; i++)
 244       histo[(int)((val[s][i] - min)/binwidth + 0.5)]++;
 245     for(i=0; i<nbin; i++)
 246       fprintf(fp," %g  %g\n",min+i*binwidth,(double)histo[i]/(n*binwidth));
 247     if (s < nset-1)
 248       fprintf(fp,"&\n");
 249   }
 250   ffclose(fp);
 251 }
 252
 253 static int real_comp(const void *a,const void *b)
 254 {
 255   real dif = *(real *)a - *(real *)b;
 256
 257   if (dif < 0)
 258     return -1;
 259   else if (dif > 0)
 260     return 1;
 261   else
 262     return 0;
 263 }
 264
 265 static void average(const char *avfile,int avbar_opt,
 266                     int n, int nset,real **val,real *t)
 267 {
 268   FILE   *fp;
 269   int    i,s,edge=0;
 270   double av,var,err;
 271   real   *tmp=NULL;
 272
 273   fp = ffopen(avfile,"w");
 274   if ((avbar_opt == avbarERROR) && (nset == 1))
 275     avbar_opt = avbarNONE;
 276   if (avbar_opt != avbarNONE) {
 277     if (avbar_opt == avbar90) {
 278       snew(tmp,nset);
 279       fprintf(fp,"@TYPE xydydy\n");
 280       edge = (int)(nset*0.05+0.5);
 281       fprintf(stdout,"Errorbars: discarding %d points on both sides: %d%%"
 282               " interval\n",edge,(int)(100*(nset-2*edge)/nset+0.5));
 283     } else
 284       fprintf(fp,"@TYPE xydy\n");
 285   }
 286
 287   for(i=0; i<n; i++) {
 288     av = 0;
 289     for(s=0; s<nset; s++)
 290       av += val[s][i];
 291     av /= nset;
 292     fprintf(fp," %g %g",t[i],av);
 293     var = 0;
 294     if (avbar_opt != avbarNONE) {
 295       if (avbar_opt == avbar90) {
 296         for(s=0; s<nset; s++)
 297           tmp[s] = val[s][i];
 298         qsort(tmp,nset,sizeof(tmp[0]),real_comp);
 299         fprintf(fp," %g %g",tmp[nset-1-edge]-av,av-tmp[edge]);
 300       } else {
 301         for(s=0; s<nset; s++)
 302           var += sqr(val[s][i]-av);
 303         if (avbar_opt == avbarSTDDEV)
 304           err = sqrt(var/nset);
 305         else
 306           err = sqrt(var/(nset*(nset-1)));
 307         fprintf(fp," %g",err);
 308       }
 309     }
 310     fprintf(fp,"\n");
 311   }
 312   ffclose(fp);
 313
 314   if (avbar_opt == avbar90)
 315     sfree(tmp);
 316 }
 317
 318 static real anal_ee_inf(real *parm,real T)
 319 {
 320   return sqrt(parm[1]*2*parm[0]/T+parm[3]*2*parm[2]/T);
 321 }
 322
 323 static real anal_ee(real *parm,real T,real t)
 324 {
 325   real e1,e2;
 326
 327   if (parm[0])
 328     e1 = exp(-t/parm[0]);
 329   else
 330     e1 = 1;
 331   if (parm[2])
 332     e2 = exp(-t/parm[2]);
 333   else
 334     e2 = 1;
 335
 336   return sqrt(parm[1]*2*parm[0]/T*((e1 - 1)*parm[0]/t + 1) +
 337               parm[3]*2*parm[2]/T*((e2 - 1)*parm[2]/t + 1));
 338 }
 339
 340 static void estimate_error(const char *eefile,int nb_min,int resol,int n,
 341                            int nset, double *av,double *sig,real **val,real dt,
 342                            gmx_bool bFitAc,gmx_bool bSingleExpFit,gmx_bool bAllowNegLTCorr,
 343                            const output_env_t oenv)
 344 {
 345     FILE   *fp;
 346     int    bs,prev_bs,nbs,nb;
 347     real   spacing,nbr;
 348     int    s,i,j;
 349     double blav,var;
 350     char   **leg;
 351     real   *tbs,*ybs,rtmp,dens,*fitsig,twooe,tau1_est,tau_sig;
 352     real   fitparm[4];
 353     real   ee,a,tau1,tau2;
 354
 355     if (n < 4)
 356     {
 357       fprintf(stdout,"The number of points is smaller than 4, can not make an error estimate\n");
 358
 359       return;
 360     }
 361
 362     fp = xvgropen(eefile,"Error estimates",
 363                   "Block size (time)","Error estimate", oenv);
 364     if (output_env_get_print_xvgr_codes(oenv))
 365     {
 366         fprintf(fp,
 367                 "@ subtitle \"using block averaging, total time %g (%d points)\"\n",
 368                 (n-1)*dt,n);
 369     }
 370     snew(leg,2*nset);
 371     xvgr_legend(fp,2*nset,(const char**)leg,oenv);
 372     sfree(leg);
 373
 374     spacing = pow(2,1.0/resol);
 375     snew(tbs,n);
 376     snew(ybs,n);
 377     snew(fitsig,n);
 378     for(s=0; s<nset; s++)
 379     {
 380         nbs = 0;
 381         prev_bs = 0;
 382         nbr = nb_min;
 383         while (nbr <= n)
 384         {
 385             bs = n/(int)nbr;
 386             if (bs != prev_bs)
 387             {
 388                 nb = n/bs;
 389                 var = 0;
 390                 for(i=0; i<nb; i++)
 391                 {
 392                     blav=0;
 393                     for (j=0; j<bs; j++)
 394                     {
 395                         blav += val[s][bs*i+j];
 396                     }
 397                     var += sqr(av[s] - blav/bs);
 398                 }
 399                 tbs[nbs] = bs*dt;
 400                 if (sig[s] == 0)
 401                 {
 402                     ybs[nbs] = 0;
 403                 }
 404                 else
 405                 {
 406                     ybs[nbs] = var/(nb*(nb-1.0))*(n*dt)/(sig[s]*sig[s]);
 407                 }
 408                 nbs++;
 409             }
 410             nbr *= spacing;
 411             nb = (int)(nbr+0.5);
 412             prev_bs = bs;
 413         }
 414         if (sig[s] == 0)
 415         {
 416             ee   = 0;
 417             a    = 1;
 418             tau1 = 0;
 419             tau2 = 0;
 420         }
 421         else
 422         {
 423             for(i=0; i<nbs/2; i++)
 424             {
 425                 rtmp         = tbs[i];
 426                 tbs[i]       = tbs[nbs-1-i];
 427                 tbs[nbs-1-i] = rtmp;
 428                 rtmp         = ybs[i];
 429                 ybs[i]       = ybs[nbs-1-i];
 430                 ybs[nbs-1-i] = rtmp;
 431             }
 432             /* The initial slope of the normalized ybs^2 is 1.
 433              * For a single exponential autocorrelation: ybs(tau1) = 2/e tau1
 434              * From this we take our initial guess for tau1.
 435              */
 436             twooe = 2/exp(1);
 437             i = -1;
 438             do
 439             {
 440                 i++;
 441                 tau1_est = tbs[i];
 442             } while (i < nbs - 1 &&
 443                      (ybs[i] > ybs[i+1] || ybs[i] > twooe*tau1_est));
 444
 445             if (ybs[0] > ybs[1])
 446             {
 447                 fprintf(stdout,"Data set %d has strange time correlations:\n"
 448                         "the std. error using single points is larger than that of blocks of 2 points\n"
 449                         "The error estimate might be inaccurate, check the fit\n",
 450                         s+1);
 451                 /* Use the total time as tau for the fitting weights */
 452                 tau_sig = (n - 1)*dt;
 453             }
 454             else
 455             {
 456                 tau_sig = tau1_est;
 457             }
 458
 459             if (debug)
 460             {
 461                 fprintf(debug,"set %d tau1 estimate %f\n",s+1,tau1_est);
 462             }
 463
 464             /* Generate more or less appropriate sigma's,
 465              * also taking the density of points into account.
 466              */
 467             for(i=0; i<nbs; i++)
 468             {
 469                 if (i == 0)
 470                 {
 471                     dens = tbs[1]/tbs[0] - 1;
 472                 }
 473                 else if (i == nbs-1)
 474                 {
 475                     dens = tbs[nbs-1]/tbs[nbs-2] - 1;
 476                 }
 477                 else
 478                 {
 479                     dens = 0.5*(tbs[i+1]/tbs[i-1] - 1);
 480                 }
 481                 fitsig[i] = sqrt((tau_sig + tbs[i])/dens);
 482             }
 483
 484             if (!bSingleExpFit)
 485             {
 486                 fitparm[0] = tau1_est;
 487                 fitparm[1] = 0.95;
 488                 /* We set the initial guess for tau2
 489                  * to halfway between tau1_est and the total time (on log scale).
 490                  */
 491                 fitparm[2] = sqrt(tau1_est*(n-1)*dt);
 492                 do_lmfit(nbs,ybs,fitsig,0,tbs,0,dt*n,oenv,
 493                          bDebugMode(),effnERREST,fitparm,0);
 494                 fitparm[3] = 1-fitparm[1];
 495             }
 496             if (bSingleExpFit || fitparm[0]<0 || fitparm[2]<0 || fitparm[1]<0
 497                 || (fitparm[1]>1 && !bAllowNegLTCorr) || fitparm[2]>(n-1)*dt)
 498             {
 499                 if (!bSingleExpFit)
 500                 {
 501                     if (fitparm[2] > (n-1)*dt)
 502                     {
 503                         fprintf(stdout,
 504                                 "Warning: tau2 is longer than the length of the data (%g)\n"
 505                                 "         the statistics might be bad\n",
 506                                 (n-1)*dt);
 507                     }
 508                     else
 509                     {
 510                         fprintf(stdout,"a fitted parameter is negative\n");
 511                     }
 512                     fprintf(stdout,"invalid fit:  e.e. %g  a %g  tau1 %g  tau2 %g\n",
 513                             sig[s]*anal_ee_inf(fitparm,n*dt),
 514                             fitparm[1],fitparm[0],fitparm[2]);
 515                     /* Do a fit with tau2 fixed at the total time.
 516                      * One could also choose any other large value for tau2.
 517                      */
 518                     fitparm[0] = tau1_est;
 519                     fitparm[1] = 0.95;
 520                     fitparm[2] = (n-1)*dt;
 521                     fprintf(stderr,"Will fix tau2 at the total time: %g\n",fitparm[2]);
 522                     do_lmfit(nbs,ybs,fitsig,0,tbs,0,dt*n,oenv,bDebugMode(),
 523                              effnERREST,fitparm,4);
 524                     fitparm[3] = 1-fitparm[1];
 525                 }
 526                 if (bSingleExpFit || fitparm[0]<0 || fitparm[1]<0
 527                     || (fitparm[1]>1 && !bAllowNegLTCorr))
 528                 {
 529                     if (!bSingleExpFit) {
 530                         fprintf(stdout,"a fitted parameter is negative\n");
 531                         fprintf(stdout,"invalid fit:  e.e. %g  a %g  tau1 %g  tau2 %g\n",
 532                                 sig[s]*anal_ee_inf(fitparm,n*dt),
 533                                 fitparm[1],fitparm[0],fitparm[2]);
 534                     }
 535                     /* Do a single exponential fit */
 536                     fprintf(stderr,"Will use a single exponential fit for set %d\n",s+1);
 537                     fitparm[0] = tau1_est;
 538                     fitparm[1] = 1.0;
 539                     fitparm[2] = 0.0;
 540                     do_lmfit(nbs,ybs,fitsig,0,tbs,0,dt*n,oenv,bDebugMode(),
 541                              effnERREST,fitparm,6);
 542                     fitparm[3] = 1-fitparm[1];
 543                 }
 544             }
 545             ee   = sig[s]*anal_ee_inf(fitparm,n*dt);
 546             a    = fitparm[1];
 547             tau1 = fitparm[0];
 548             tau2 = fitparm[2];
 549         }
 550         fprintf(stdout,"Set %3d:  err.est. %g  a %g  tau1 %g  tau2 %g\n",
 551                 s+1,ee,a,tau1,tau2);
 552         fprintf(fp,"@ legend string %d \"av %f\"\n",2*s,av[s]);
 553         fprintf(fp,"@ legend string %d \"ee %6g\"\n",
 554                 2*s+1,sig[s]*anal_ee_inf(fitparm,n*dt));
 555         for(i=0; i<nbs; i++)
 556         {
 557             fprintf(fp,"%g %g %g\n",tbs[i],sig[s]*sqrt(ybs[i]/(n*dt)),
 558                     sig[s]*sqrt(fit_function(effnERREST,fitparm,tbs[i])/(n*dt)));
 559         }
 560
 561         if (bFitAc)
 562         {
 563             int fitlen;
 564             real *ac,acint,ac_fit[4];
 565
 566             snew(ac,n);
 567             for(i=0; i<n; i++) {
 568                 ac[i] = val[s][i] - av[s];
 569                 if (i > 0)
 570                     fitsig[i] = sqrt(i);
 571                 else
 572                     fitsig[i] = 1;
 573             }
 574             low_do_autocorr(NULL,oenv,NULL,n,1,-1,&ac,
 575                             dt,eacNormal,1,FALSE,TRUE,
 576                             FALSE,0,0,effnNONE,0);
 577
 578             fitlen = n/nb_min;
 579
 580             /* Integrate ACF only up to fitlen/2 to avoid integrating noise */
 581             acint = 0.5*ac[0];
 582             for(i=1; i<=fitlen/2; i++)
 583             {
 584                 acint += ac[i];
 585             }
 586             acint *= dt;
 587
 588             /* Generate more or less appropriate sigma's */
 589             for(i=0; i<=fitlen; i++)
 590             {
 591                 fitsig[i] = sqrt(acint + dt*i);
 592             }
 593
 594             ac_fit[0] = 0.5*acint;
 595             ac_fit[1] = 0.95;
 596             ac_fit[2] = 10*acint;
 597             do_lmfit(n/nb_min,ac,fitsig,dt,0,0,fitlen*dt,oenv,
 598                      bDebugMode(),effnEXP3,ac_fit,0);
 599             ac_fit[3] = 1 - ac_fit[1];
 600
 601             fprintf(stdout,"Set %3d:  ac erest %g  a %g  tau1 %g  tau2 %g\n",
 602                     s+1,sig[s]*anal_ee_inf(ac_fit,n*dt),
 603                     ac_fit[1],ac_fit[0],ac_fit[2]);
 604
 605             fprintf(fp,"&\n");
 606             for(i=0; i<nbs; i++)
 607             {
 608                 fprintf(fp,"%g %g\n",tbs[i],
 609                         sig[s]*sqrt(fit_function(effnERREST,ac_fit,tbs[i]))/(n*dt));
 610             }
 611
 612             sfree(ac);
 613         }
 614         if (s < nset-1)
 615         {
 616             fprintf(fp,"&\n");
 617         }
 618     }
 619     sfree(fitsig);
 620     sfree(ybs);
 621     sfree(tbs);
 622     ffclose(fp);
 623 }
 624
 625 static void luzar_correl(int nn,real *time,int nset,real **val,real temp,
 626                          gmx_bool bError,real fit_start,real smooth_tail_start,
 627                          const output_env_t oenv)
 628 {
 629   const real tol = 1e-8;
 630   real *kt;
 631   real weight,d2;
 632   int  j;
 633
 634   please_cite(stdout,"Spoel2006b");
 635
 636   /* Compute negative derivative k(t) = -dc(t)/dt */
 637   if (!bError) {
 638     snew(kt,nn);
 639     compute_derivative(nn,time,val[0],kt);
 640     for(j=0; (j<nn); j++)
 641       kt[j] = -kt[j];
 642     if (debug) {
 643       d2 = 0;
 644       for(j=0; (j<nn); j++)
 645         d2 += sqr(kt[j] - val[3][j]);
 646       fprintf(debug,"RMS difference in derivatives is %g\n",sqrt(d2/nn));
 647     }
 648     analyse_corr(nn,time,val[0],val[2],kt,NULL,NULL,NULL,fit_start,
 649                  temp,smooth_tail_start,oenv);
 650     sfree(kt);
 651   }
 652   else if (nset == 6) {
 653     analyse_corr(nn,time,val[0],val[2],val[4],
 654                  val[1],val[3],val[5],fit_start,temp,smooth_tail_start,oenv);
 655   }
 656   else {
 657     printf("Inconsistent input. I need c(t) sigma_c(t) n(t) sigma_n(t) K(t) sigma_K(t)\n");
 658     printf("Not doing anything. Sorry.\n");
 659   }
 660 }
 661
 662 static void filter(real flen,int n,int nset,real **val,real dt,
 663                    const output_env_t oenv)
 664 {
 665   int    f,s,i,j;
 666   double *filt,sum,vf,fluc,fluctot;
 667
 668   f = (int)(flen/(2*dt));
 669   snew(filt,f+1);
 670   filt[0] = 1;
 671   sum = 1;
 672   for(i=1; i<=f; i++) {
 673     filt[i] = cos(M_PI*dt*i/flen);
 674     sum += 2*filt[i];
 675   }
 676   for(i=0; i<=f; i++)
 677     filt[i] /= sum;
 678   fprintf(stdout,"Will calculate the fluctuation over %d points\n",n-2*f);
 679   fprintf(stdout,"  using a filter of length %g of %d points\n",flen,2*f+1);
 680   fluctot = 0;
 681   for(s=0; s<nset; s++) {
 682     fluc = 0;
 683     for(i=f; i<n-f; i++) {
 684       vf = filt[0]*val[s][i];
 685       for(j=1; j<=f; j++)
 686         vf += filt[j]*(val[s][i-f]+val[s][i+f]);
 687       fluc += sqr(val[s][i] - vf);
 688     }
 689     fluc /= n - 2*f;
 690     fluctot += fluc;
 691     fprintf(stdout,"Set %3d filtered fluctuation: %12.6e\n",s+1,sqrt(fluc));
 692   }
 693   fprintf(stdout,"Overall filtered fluctuation: %12.6e\n",sqrt(fluctot/nset));
 694   fprintf(stdout,"\n");
 695
 696   sfree(filt);
 697 }
 698
 699 static void do_fit(FILE *out,int n,gmx_bool bYdy,
 700                    int ny,real *x0,real **val,
 701                    int npargs,t_pargs *ppa,const output_env_t oenv)
 702 {
 703   real *c1=NULL,*sig=NULL,*fitparm;
 704   real tendfit,tbeginfit;
 705   int  i,efitfn,nparm;
 706
 707   efitfn = get_acffitfn();
 708   nparm  = nfp_ffn[efitfn];
 709   fprintf(out,"Will fit to the following function:\n");
 710   fprintf(out,"%s\n",longs_ffn[efitfn]);
 711   c1 = val[n];
 712   if (bYdy) {
 713     c1  = val[n];
 714     sig = val[n+1];
 715     fprintf(out,"Using two columns as y and sigma values\n");
 716   } else {
 717     snew(sig,ny);
 718   }
 719   if (opt2parg_bSet("-beginfit",npargs,ppa)) {
 720     tbeginfit = opt2parg_real("-beginfit",npargs,ppa);
 721   } else {
 722     tbeginfit = x0[0];
 723   }
 724   if (opt2parg_bSet("-endfit",npargs,ppa)) {
 725     tendfit   = opt2parg_real("-endfit",npargs,ppa);
 726   } else {
 727     tendfit   = x0[ny-1];
 728   }
 729
 730   snew(fitparm,nparm);
 731   switch(efitfn) {
 732   case effnEXP1:
 733     fitparm[0] = 0.5;
 734     break;
 735   case effnEXP2:
 736     fitparm[0] = 0.5;
 737     fitparm[1] = c1[0];
 738     break;
 739   case effnEXP3:
 740     fitparm[0] = 1.0;
 741     fitparm[1] = 0.5*c1[0];
 742     fitparm[2] = 10.0;
 743     break;
 744   case effnEXP5:
 745     fitparm[0] = fitparm[2] = 0.5*c1[0];
 746     fitparm[1] = 10;
 747     fitparm[3] = 40;
 748     fitparm[4] = 0;
 749     break;
 750   case effnEXP7:
 751     fitparm[0] = fitparm[2] = fitparm[4] = 0.33*c1[0];
 752     fitparm[1] = 1;
 753     fitparm[3] = 10;
 754     fitparm[5] = 100;
 755     fitparm[6] = 0;
 756     break;
 757   case effnEXP9:
 758     fitparm[0] = fitparm[2] = fitparm[4] = fitparm[6] = 0.25*c1[0];
 759     fitparm[1] = 0.1;
 760     fitparm[3] = 1;
 761     fitparm[5] = 10;
 762     fitparm[7] = 100;
 763     fitparm[8] = 0;
 764     break;
 765   default:
 766     fprintf(out,"Warning: don't know how to initialize the parameters\n");
 767     for(i=0; (i<nparm); i++)
 768       fitparm[i] = 1;
 769   }
 770   fprintf(out,"Starting parameters:\n");
 771   for(i=0; (i<nparm); i++)
 772     fprintf(out,"a%-2d = %12.5e\n",i+1,fitparm[i]);
 773   if (do_lmfit(ny,c1,sig,0,x0,tbeginfit,tendfit,
 774                oenv,bDebugMode(),efitfn,fitparm,0)) {
 775     for(i=0; (i<nparm); i++)
 776       fprintf(out,"a%-2d = %12.5e\n",i+1,fitparm[i]);
 777   }
 778   else {
 779     fprintf(out,"No solution was found\n");
 780   }
 781 }
 782
 783 static void do_ballistic(const char *balFile, int nData,
 784                          real *t, real **val, int nSet,
 785                          real balTime, int nBalExp,
 786                          gmx_bool bDerivative,
 787                          const output_env_t oenv)
 788 {
 789   double **ctd=NULL, *td=NULL;
 790   t_gemParams *GP = init_gemParams(0, 0, t, nData, 0, 0, 0, balTime, nBalExp, bDerivative);
 791   static char *leg[] = {"Ac'(t)"};
 792   FILE *fp;
 793   int i, set;
 794
 795   if (GP->ballistic/GP->tDelta >= GP->nExpFit*2+1)
 796   {
 797       snew(ctd, nSet);
 798       snew(td,  nData);
 799
 800       fp = xvgropen(balFile, "Hydrogen Bond Autocorrelation","Time (ps)","C'(t)", oenv);
 801       xvgr_legend(fp,asize(leg),(const char**)leg,oenv);
 802
 803       for (set=0; set<nSet; set++)
 804       {
 805           snew(ctd[set], nData);
 806           for (i=0; i<nData; i++) {
 807               ctd[set][i] = (double)val[set][i];
 808               if (set==0)
 809                   td[i] = (double)t[i];
 810           }
 811
 812           takeAwayBallistic(ctd[set], td, nData, GP->ballistic, GP->nExpFit, GP->bDt);
 813       }
 814
 815       for (i=0; i<nData; i++)
 816       {
 817           fprintf(fp, "  %g",t[i]);
 818           for (set=0; set<nSet; set++)
 819           {
 820               fprintf(fp, "  %g", ctd[set][i]);
 821           }
 822           fprintf(fp, "\n");
 823       }
 824
 825
 826       for (set=0; set<nSet; set++)
 827           sfree(ctd[set]);
 828       sfree(ctd);
 829       sfree(td);
 830   }
 831   else
 832       printf("Number of data points is less than the number of parameters to fit\n."
 833              "The system is underdetermined, hence no ballistic term can be found.\n\n");
 834 }
 835
 836 static void do_geminate(const char *gemFile, int nData,
 837                         real *t, real **val, int nSet,
 838                         const real D, const real rcut, const real balTime,
 839                         const int nFitPoints, const real begFit, const real endFit,
 840                         const output_env_t oenv)
 841 {
 842     double **ctd=NULL, **ctdGem=NULL, *td=NULL;
 843     t_gemParams *GP = init_gemParams(rcut, D, t, nData, nFitPoints,
 844                                      begFit, endFit, balTime, 1, FALSE);
 845     const char *leg[] = {"Ac\\sgem\\N(t)"};
 846     FILE *fp;
 847     int i, set;
 848
 849     snew(ctd,    nSet);
 850     snew(ctdGem, nSet);
 851     snew(td,  nData);
 852
 853     fp = xvgropen(gemFile, "Hydrogen Bond Autocorrelation","Time (ps)","C'(t)", oenv);
 854     xvgr_legend(fp,asize(leg),leg,oenv);
 855
 856     for (set=0; set<nSet; set++)
 857     {
 858         snew(ctd[set],    nData);
 859         snew(ctdGem[set], nData);
 860         for (i=0; i<nData; i++) {
 861             ctd[set][i] = (double)val[set][i];
 862             if (set==0)
 863                 td[i] = (double)t[i];
 864         }
 865         fitGemRecomb(ctd[set], td, &(ctd[set]), nData, GP);
 866     }
 867
 868     for (i=0; i<nData; i++)
 869         {
 870         fprintf(fp, "  %g",t[i]);
 871         for (set=0; set<nSet; set++)
 872             {
 873             fprintf(fp, "  %g", ctdGem[set][i]);
 874             }
 875         fprintf(fp, "\n");
 876         }
 877
 878     for (set=0; set<nSet; set++)
 879     {
 880         sfree(ctd[set]);
 881         sfree(ctdGem[set]);
 882     }
 883     sfree(ctd);
 884     sfree(ctdGem);
 885     sfree(td);
 886 }
 887
 888 int gmx_analyze(int argc,char *argv[])
 889 {
 890   static const char *desc[] = {
 891     "[TT]g_analyze[tt] reads an ascii file and analyzes data sets.",
 892     "A line in the input file may start with a time",
 893     "(see option [TT]-time[tt]) and any number of y values may follow.",
 894     "Multiple sets can also be",
 895     "read when they are separated by & (option [TT]-n[tt]),",
 896     "in this case only one y value is read from each line.",
 897     "All lines starting with # and @ are skipped.",
 898     "All analyses can also be done for the derivative of a set",
 899     "(option [TT]-d[tt]).[PAR]",
 900
 901     "All options, except for [TT]-av[tt] and [TT]-power[tt] assume that the",
 902     "points are equidistant in time.[PAR]",
 903
 904     "[TT]g_analyze[tt] always shows the average and standard deviation of each",
 905     "set. For each set it also shows the relative deviation of the third",
 906     "and fourth cumulant from those of a Gaussian distribution with the same",
 907     "standard deviation.[PAR]",
 908
 909     "Option [TT]-ac[tt] produces the autocorrelation function(s).[PAR]",
 910
 911     "Option [TT]-cc[tt] plots the resemblance of set i with a cosine of",
 912     "i/2 periods. The formula is:[BR]"
 913     "2 (int0-T y(t) cos(i pi t) dt)^2 / int0-T y(t) y(t) dt[BR]",
 914     "This is useful for principal components obtained from covariance",
 915     "analysis, since the principal components of random diffusion are",
 916     "pure cosines.[PAR]",
 917
 918     "Option [TT]-msd[tt] produces the mean square displacement(s).[PAR]",
 919
 920     "Option [TT]-dist[tt] produces distribution plot(s).[PAR]",
 921
 922     "Option [TT]-av[tt] produces the average over the sets.",
 923     "Error bars can be added with the option [TT]-errbar[tt].",
 924     "The errorbars can represent the standard deviation, the error",
 925     "(assuming the points are independent) or the interval containing",
 926     "90% of the points, by discarding 5% of the points at the top and",
 927     "the bottom.[PAR]",
 928
 929     "Option [TT]-ee[tt] produces error estimates using block averaging.",
 930     "A set is divided in a number of blocks and averages are calculated for",
 931     "each block. The error for the total average is calculated from",
 932     "the variance between averages of the m blocks B_i as follows:",
 933     "error^2 = Sum (B_i - <B>)^2 / (m*(m-1)).",
 934     "These errors are plotted as a function of the block size.",
 935     "Also an analytical block average curve is plotted, assuming",
 936     "that the autocorrelation is a sum of two exponentials.",
 937     "The analytical curve for the block average is:[BR]",
 938     "f(t) = sigma sqrt(2/T (  a   (tau1 ((exp(-t/tau1) - 1) tau1/t + 1)) +[BR]",
 939     "                       (1-a) (tau2 ((exp(-t/tau2) - 1) tau2/t + 1)))),[BR]"
 940     "where T is the total time.",
 941     "a, tau1 and tau2 are obtained by fitting f^2(t) to error^2.",
 942     "When the actual block average is very close to the analytical curve,",
 943     "the error is sigma*sqrt(2/T (a tau1 + (1-a) tau2)).",
 944     "The complete derivation is given in",
 945     "B. Hess, J. Chem. Phys. 116:209-217, 2002.[PAR]",
 946
 947     "Option [TT]-bal[tt] finds and subtracts the ultrafast \"ballistic\"",
 948     "component from a hydrogen bond autocorrelation function by the fitting",
 949     "of a sum of exponentials, as described in e.g.",
 950     "O. Markovitch, J. Chem. Phys. 129:084505, 2008. The fastest term",
 951     "is the one with the most negative coefficient in the exponential,",
 952     "or with [TT]-d[tt], the one with most negative time derivative at time 0.",
 953     "[TT]-nbalexp[tt] sets the number of exponentials to fit.[PAR]",
 954
 955     "Option [TT]-gem[tt] fits bimolecular rate constants ka and kb",
 956     "(and optionally kD) to the hydrogen bond autocorrelation function",
 957     "according to the reversible geminate recombination model. Removal of",
 958     "the ballistic component first is strongly adviced. The model is presented in",
 959     "O. Markovitch, J. Chem. Phys. 129:084505, 2008.[PAR]",
 960
 961     "Option [TT]-filter[tt] prints the RMS high-frequency fluctuation",
 962     "of each set and over all sets with respect to a filtered average.",
 963     "The filter is proportional to cos(pi t/len) where t goes from -len/2",
 964     "to len/2. len is supplied with the option [TT]-filter[tt].",
 965     "This filter reduces oscillations with period len/2 and len by a factor",
 966     "of 0.79 and 0.33 respectively.[PAR]",
 967
 968     "Option [TT]-g[tt] fits the data to the function given with option",
 969     "[TT]-fitfn[tt].[PAR]",
 970
 971     "Option [TT]-power[tt] fits the data to b t^a, which is accomplished",
 972     "by fitting to a t + b on log-log scale. All points after the first",
 973     "zero or negative value are ignored.[PAR]"
 974
 975     "Option [TT]-luzar[tt] performs a Luzar & Chandler kinetics analysis",
 976     "on output from [TT]g_hbond[tt]. The input file can be taken directly",
 977     "from [TT]g_hbond -ac[tt], and then the same result should be produced."
 978   };
 979   static real tb=-1,te=-1,frac=0.5,filtlen=0,binwidth=0.1,aver_start=0;
 980   static gmx_bool bHaveT=TRUE,bDer=FALSE,bSubAv=TRUE,bAverCorr=FALSE,bXYdy=FALSE;
 981   static gmx_bool bEESEF=FALSE,bEENLC=FALSE,bEeFitAc=FALSE,bPower=FALSE;
 982   static gmx_bool bIntegrate=FALSE,bRegression=FALSE,bLuzar=FALSE,bLuzarError=FALSE;
 983   static int  nsets_in=1,d=1,nb_min=4,resol=10, nBalExp=4, nFitPoints=100;
 984   static real temp=298.15,fit_start=1, fit_end=60, smooth_tail_start=-1, balTime=0.2, diffusion=5e-5,rcut=0.35;
 985
 986   /* must correspond to enum avbar* declared at beginning of file */
 987   static const char *avbar_opt[avbarNR+1] = {
 988     NULL, "none", "stddev", "error", "90", NULL
 989   };
 990
 991   t_pargs pa[] = {
 992     { "-time",    FALSE, etBOOL, {&bHaveT},
 993       "Expect a time in the input" },
 994     { "-b",       FALSE, etREAL, {&tb},
 995       "First time to read from set" },
 996     { "-e",       FALSE, etREAL, {&te},
 997       "Last time to read from set" },
 998     { "-n",       FALSE, etINT, {&nsets_in},
 999       "Read # sets separated by &" },
1000     { "-d",       FALSE, etBOOL, {&bDer},
1001         "Use the derivative" },
1002     { "-dp",      FALSE, etINT, {&d},
1003       "HIDDENThe derivative is the difference over # points" },
1004     { "-bw",      FALSE, etREAL, {&binwidth},
1005       "Binwidth for the distribution" },
1006     { "-errbar",  FALSE, etENUM, {avbar_opt},
1007       "Error bars for [TT]-av[tt]" },
1008     { "-integrate",FALSE,etBOOL, {&bIntegrate},
1009       "Integrate data function(s) numerically using trapezium rule" },
1010     { "-aver_start",FALSE, etREAL, {&aver_start},
1011       "Start averaging the integral from here" },
1012     { "-xydy",    FALSE, etBOOL, {&bXYdy},
1013       "Interpret second data set as error in the y values for integrating" },
1014     { "-regression",FALSE,etBOOL,{&bRegression},
1015       "Perform a linear regression analysis on the data. If [TT]-xydy[tt] is set a second set will be interpreted as the error bar in the Y value. Otherwise, if multiple data sets are present a multilinear regression will be performed yielding the constant A that minimize chi^2 = (y - A0 x0 - A1 x1 - ... - AN xN)^2 where now Y is the first data set in the input file and xi the others. Do read the information at the option [TT]-time[tt]." },
1016     { "-luzar",   FALSE, etBOOL, {&bLuzar},
1017       "Do a Luzar and Chandler analysis on a correlation function and related as produced by [TT]g_hbond[tt]. When in addition the [TT]-xydy[tt] flag is given the second and fourth column will be interpreted as errors in c(t) and n(t)." },
1018     { "-temp",    FALSE, etREAL, {&temp},
1019       "Temperature for the Luzar hydrogen bonding kinetics analysis" },
1020     { "-fitstart", FALSE, etREAL, {&fit_start},
1021       "Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation" },
1022     { "-fitend", FALSE, etREAL, {&fit_end},
1023       "Time (ps) where to stop fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation. Only with [TT]-gem[tt]" },
1024     { "-smooth",FALSE, etREAL, {&smooth_tail_start},
1025       "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/tau)" },
1026     { "-nbmin",   FALSE, etINT, {&nb_min},
1027       "HIDDENMinimum number of blocks for block averaging" },
1028     { "-resol", FALSE, etINT, {&resol},
1029       "HIDDENResolution for the block averaging, block size increases with"
1030     " a factor 2^(1/#)" },
1031     { "-eeexpfit", FALSE, etBOOL, {&bEESEF},
1032       "HIDDENAlways use a single exponential fit for the error estimate" },
1033     { "-eenlc", FALSE, etBOOL, {&bEENLC},
1034       "HIDDENAllow a negative long-time correlation" },
1035     { "-eefitac", FALSE, etBOOL, {&bEeFitAc},
1036       "HIDDENAlso plot analytical block average using a autocorrelation fit" },
1037     { "-filter",  FALSE, etREAL, {&filtlen},
1038       "Print the high-frequency fluctuation after filtering with a cosine filter of length #" },
1039     { "-power", FALSE, etBOOL, {&bPower},
1040       "Fit data to: b t^a" },
1041     { "-subav", FALSE, etBOOL, {&bSubAv},
1042       "Subtract the average before autocorrelating" },
1043     { "-oneacf", FALSE, etBOOL, {&bAverCorr},
1044       "Calculate one ACF over all sets" },
1045     { "-nbalexp", FALSE, etINT, {&nBalExp},
1046       "HIDDENNumber of exponentials to fit to the ultrafast component" },
1047     { "-baltime", FALSE, etREAL, {&balTime},
1048       "HIDDENTime up to which the ballistic component will be fitted" },
1049 /*     { "-gemnp", FALSE, etINT, {&nFitPoints}, */
1050 /*       "HIDDENNumber of data points taken from the ACF to use for fitting to rev. gem. recomb. model."}, */
1051 /*     { "-rcut", FALSE, etREAL, {&rcut}, */
1052 /*       "Cut-off for hydrogen bonds in geminate algorithms" }, */
1053 /*     { "-gemtype", FALSE, etENUM, {gemType}, */
1054 /*       "What type of gminate recombination to use"}, */
1055 /*     { "-D", FALSE, etREAL, {&diffusion}, */
1056 /*       "The self diffusion coefficient which is used for the reversible geminate recombination model."} */
1057   };
1058 #define NPA asize(pa)
1059
1060   FILE     *out,*out_fit;
1061   int      n,nlast,s,nset,i,j=0;
1062   real     **val,*t,dt,tot,error;
1063   double   *av,*sig,cum1,cum2,cum3,cum4,db;
1064   const char     *acfile,*msdfile,*ccfile,*distfile,*avfile,*eefile,*balfile,*gemfile,*fitfile;
1065   output_env_t oenv;
1066
1067   t_filenm fnm[] = {
1068     { efXVG, "-f",    "graph",    ffREAD   },
1069     { efXVG, "-ac",   "autocorr", ffOPTWR  },
1070     { efXVG, "-msd",  "msd",      ffOPTWR  },
1071     { efXVG, "-cc",   "coscont",  ffOPTWR  },
1072     { efXVG, "-dist", "distr",    ffOPTWR  },
1073     { efXVG, "-av",   "average",  ffOPTWR  },
1074     { efXVG, "-ee",   "errest",   ffOPTWR  },
1075     { efXVG, "-bal",  "ballisitc",ffOPTWR  },
1076 /*     { efXVG, "-gem",  "geminate", ffOPTWR  }, */
1077     { efLOG, "-g",    "fitlog",   ffOPTWR  }
1078   };
1079 #define NFILE asize(fnm)
1080
1081   int     npargs;
1082   t_pargs *ppa;
1083
1084   npargs = asize(pa);
1085   ppa    = add_acf_pargs(&npargs,pa);
1086
1087   CopyRight(stderr,argv[0]);
1088   parse_common_args(&argc,argv,PCA_CAN_VIEW,
1089                     NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1090
1091   acfile   = opt2fn_null("-ac",NFILE,fnm);
1092   msdfile  = opt2fn_null("-msd",NFILE,fnm);
1093   ccfile   = opt2fn_null("-cc",NFILE,fnm);
1094   distfile = opt2fn_null("-dist",NFILE,fnm);
1095   avfile   = opt2fn_null("-av",NFILE,fnm);
1096   eefile   = opt2fn_null("-ee",NFILE,fnm);
1097   balfile  = opt2fn_null("-bal",NFILE,fnm);
1098 /*   gemfile  = opt2fn_null("-gem",NFILE,fnm); */
1099     /* When doing autocorrelation we don't want a fitlog for fitting
1100      * the function itself (not the acf) when the user did not ask for it.
1101      */
1102     if (opt2parg_bSet("-fitfn",npargs,ppa) && acfile == NULL)
1103     {
1104         fitfile  = opt2fn("-g",NFILE,fnm);
1105     }
1106     else
1107     {
1108         fitfile  = opt2fn_null("-g",NFILE,fnm);
1109     }
1110
1111     val = read_xvg_time(opt2fn("-f",NFILE,fnm),bHaveT,
1112                         opt2parg_bSet("-b",npargs,ppa),tb,
1113                         opt2parg_bSet("-e",npargs,ppa),te,
1114                         nsets_in,&nset,&n,&dt,&t);
1115     printf("Read %d sets of %d points, dt = %g\n\n",nset,n,dt);
1116
1117     if (bDer)
1118     {
1119         printf("Calculating the derivative as (f[i+%d]-f[i])/(%d*dt)\n\n",
1120                d,d);
1121         n -= d;
1122         for(s=0; s<nset; s++)
1123         {
1124             for(i=0; (i<n); i++)
1125             {
1126                 val[s][i] = (val[s][i+d]-val[s][i])/(d*dt);
1127             }
1128         }
1129     }
1130
1131     if (bIntegrate)
1132     {
1133         real sum,stddev;
1134
1135         printf("Calculating the integral using the trapezium rule\n");
1136
1137         if (bXYdy)
1138         {
1139             sum = evaluate_integral(n,t,val[0],val[1],aver_start,&stddev);
1140             printf("Integral %10.3f +/- %10.5f\n",sum,stddev);
1141         }
1142         else
1143         {
1144             for(s=0; s<nset; s++)
1145             {
1146                 sum = evaluate_integral(n,t,val[s],NULL,aver_start,&stddev);
1147                 printf("Integral %d  %10.5f  +/- %10.5f\n",s+1,sum,stddev);
1148             }
1149         }
1150     }
1151
1152     if (fitfile != NULL)
1153     {
1154         out_fit = ffopen(fitfile,"w");
1155         if (bXYdy && nset >= 2)
1156         {
1157             do_fit(out_fit,0,TRUE,n,t,val,npargs,ppa,oenv);
1158         }
1159         else
1160         {
1161             for(s=0; s<nset; s++)
1162             {
1163                 do_fit(out_fit,s,FALSE,n,t,val,npargs,ppa,oenv);
1164             }
1165         }
1166         ffclose(out_fit);
1167     }
1168
1169   printf("                                      std. dev.    relative deviation of\n");
1170   printf("                       standard       ---------   cumulants from those of\n");
1171   printf("set      average       deviation      sqrt(n-1)   a Gaussian distribition\n");
1172   printf("                                                      cum. 3   cum. 4\n");
1173   snew(av,nset);
1174   snew(sig,nset);
1175   for(s=0; (s<nset); s++) {
1176     cum1 = 0;
1177     cum2 = 0;
1178     cum3 = 0;
1179     cum4 = 0;
1180     for(i=0; (i<n); i++)
1181       cum1 += val[s][i];
1182     cum1 /= n;
1183     for(i=0; (i<n); i++) {
1184       db = val[s][i]-cum1;
1185       cum2 += db*db;
1186       cum3 += db*db*db;
1187       cum4 += db*db*db*db;
1188     }
1189     cum2  /= n;
1190     cum3  /= n;
1191     cum4  /= n;
1192     av[s]  = cum1;
1193     sig[s] = sqrt(cum2);
1194     if (n > 1)
1195       error = sqrt(cum2/(n-1));
1196     else
1197       error = 0;
1198     printf("SS%d  %13.6e   %12.6e   %12.6e      %6.3f   %6.3f\n",
1199            s+1,av[s],sig[s],error,
1200            sig[s] ? cum3/(sig[s]*sig[s]*sig[s]*sqrt(8/M_PI)) : 0,
1201            sig[s] ? cum4/(sig[s]*sig[s]*sig[s]*sig[s]*3)-1 : 0);
1202   }
1203   printf("\n");
1204
1205   if (filtlen)
1206     filter(filtlen,n,nset,val,dt,oenv);
1207
1208   if (msdfile) {
1209     out=xvgropen(msdfile,"Mean square displacement",
1210                  "time","MSD (nm\\S2\\N)",oenv);
1211     nlast = (int)(n*frac);
1212     for(s=0; s<nset; s++) {
1213       for(j=0; j<=nlast; j++) {
1214         if (j % 100 == 0)
1215           fprintf(stderr,"\r%d",j);
1216         tot=0;
1217         for(i=0; i<n-j; i++)
1218           tot += sqr(val[s][i]-val[s][i+j]);
1219         tot /= (real)(n-j);
1220         fprintf(out," %g %8g\n",dt*j,tot);
1221       }
1222       if (s<nset-1)
1223         fprintf(out,"&\n");
1224     }
1225     ffclose(out);
1226     fprintf(stderr,"\r%d, time=%g\n",j-1,(j-1)*dt);
1227   }
1228   if (ccfile)
1229     plot_coscont(ccfile,n,nset,val,oenv);
1230
1231   if (distfile)
1232     histogram(distfile,binwidth,n,nset,val,oenv);
1233   if (avfile)
1234     average(avfile,nenum(avbar_opt),n,nset,val,t);
1235   if (eefile)
1236     estimate_error(eefile,nb_min,resol,n,nset,av,sig,val,dt,
1237                    bEeFitAc,bEESEF,bEENLC,oenv);
1238   if (balfile)
1239       do_ballistic(balfile,n,t,val,nset,balTime,nBalExp,bDer,oenv);
1240 /*   if (gemfile) */
1241 /*       do_geminate(gemfile,n,t,val,nset,diffusion,rcut,balTime, */
1242 /*                   nFitPoints, fit_start, fit_end, oenv); */
1243   if (bPower)
1244     power_fit(n,nset,val,t);
1245
1246     if (acfile != NULL)
1247     {
1248         if (bSubAv)
1249         {
1250             for(s=0; s<nset; s++)
1251             {
1252                 for(i=0; i<n; i++)
1253                 {
1254                     val[s][i] -= av[s];
1255                 }
1256             }
1257         }
1258         do_autocorr(acfile,oenv,"Autocorrelation",n,nset,val,dt,
1259                     eacNormal,bAverCorr);
1260     }
1261
1262   if (bRegression)
1263       regression_analysis(n,bXYdy,t,nset,val);
1264
1265   if (bLuzar)
1266     luzar_correl(n,t,nset,val,temp,bXYdy,fit_start,smooth_tail_start,oenv);
1267
1268   view_all(oenv,NFILE, fnm);
1269
1270   thanx(stderr);
1271
1272   return 0;
1273 }
1274