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