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