src/gromacs/gmxana/gmx_analyze.cpp

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