src/gromacs/gmxana/gmx_analyze.c

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