src/gromacs/correlationfunctions/autocorr.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, 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 /*! \internal \file
  39  * \brief
  40  * Implements function to compute many autocorrelation functions
  41  *
  42  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  43  * \ingroup module_correlationfunctions
  44  */
  45 #include "gmxpre.h"
  46
  47 #include "autocorr.h"
  48
  49 #include <cmath>
  50 #include <cstdio>
  51 #include <cstring>
  52
  53 #include <algorithm>
  54
  55 #include "gromacs/correlationfunctions/expfit.h"
  56 #include "gromacs/correlationfunctions/integrate.h"
  57 #include "gromacs/correlationfunctions/manyautocorrelation.h"
  58 #include "gromacs/correlationfunctions/polynomials.h"
  59 #include "gromacs/fileio/xvgr.h"
  60 #include "gromacs/math/functions.h"
  61 #include "gromacs/math/vec.h"
  62 #include "gromacs/utility/arraysize.h"
  63 #include "gromacs/utility/fatalerror.h"
  64 #include "gromacs/utility/futil.h"
  65 #include "gromacs/utility/real.h"
  66 #include "gromacs/utility/smalloc.h"
  67 #include "gromacs/utility/strconvert.h"
  68
  69 /*! \brief Shortcut macro to select modes. */
  70 #define MODE(x) ((mode & (x)) == (x))
  71
  72 typedef struct
  73 {
  74     unsigned long mode;
  75     int           nrestart, nout, P, fitfn;
  76     gmx_bool      bFour, bNormalize;
  77     real          tbeginfit, tendfit;
  78 } t_acf;
  79
  80 /*! \brief Global variable set true if initialization routines are called. */
  81 static gmx_bool bACFinit = FALSE;
  82
  83 /*! \brief Data structure for storing command line variables. */
  84 static t_acf acf;
  85
  86 enum
  87 {
  88     enNorm,
  89     enCos,
  90     enSin
  91 };
  92
  93 /*! \brief Routine to compute ACF using FFT. */
  94 static void low_do_four_core(int nframes, real c1[], real cfour[], int nCos)
  95 {
  96     int                            i = 0;
  97     std::vector<std::vector<real>> data;
  98     data.resize(1);
  99     data[0].resize(nframes, 0);
 100     switch (nCos)
 101     {
 102         case enNorm:
 103             for (i = 0; (i < nframes); i++)
 104             {
 105                 data[0][i] = c1[i];
 106             }
 107             break;
 108         case enCos:
 109             for (i = 0; (i < nframes); i++)
 110             {
 111                 data[0][i] = cos(c1[i]);
 112             }
 113             break;
 114         case enSin:
 115             for (i = 0; (i < nframes); i++)
 116             {
 117                 data[0][i] = sin(c1[i]);
 118             }
 119             break;
 120         default: gmx_fatal(FARGS, "nCos = %d, %s %d", nCos, __FILE__, __LINE__);
 121     }
 122
 123     many_auto_correl(&data);
 124     for (i = 0; (i < nframes); i++)
 125     {
 126         cfour[i] = data[0][i];
 127     }
 128 }
 129
 130 /*! \brief Routine to comput ACF without FFT. */
 131 static void do_ac_core(int nframes, int nout, real corr[], real c1[], int nrestart, unsigned long mode)
 132 {
 133     int  j, k, j3, jk3, m, n;
 134     real ccc, cth;
 135     rvec xj, xk;
 136
 137     if (nrestart < 1)
 138     {
 139         printf("WARNING: setting number of restarts to 1\n");
 140         nrestart = 1;
 141     }
 142     if (debug)
 143     {
 144         fprintf(debug, "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n", nframes, nout, nrestart, mode);
 145     }
 146
 147     for (j = 0; (j < nout); j++)
 148     {
 149         corr[j] = 0;
 150     }
 151
 152     /* Loop over starting points. */
 153     for (j = 0; (j < nframes); j += nrestart)
 154     {
 155         j3 = DIM * j;
 156
 157         /* Loop over the correlation length for this starting point */
 158         for (k = 0; (k < nout) && (j + k < nframes); k++)
 159         {
 160             jk3 = DIM * (j + k);
 161
 162             /* Switch over possible ACF types.
 163              * It might be more efficient to put the loops inside the switch,
 164              * but this is more clear, and save development time!
 165              */
 166             if (MODE(eacNormal))
 167             {
 168                 corr[k] += c1[j] * c1[j + k];
 169             }
 170             else if (MODE(eacCos))
 171             {
 172                 /* Compute the cos (phi(t)-phi(t+dt)) */
 173                 corr[k] += std::cos(c1[j] - c1[j + k]);
 174             }
 175             else if (MODE(eacIden))
 176             {
 177                 /* Check equality (phi(t)==phi(t+dt)) */
 178                 corr[k] += (c1[j] == c1[j + k]) ? 1 : 0;
 179             }
 180             else if (MODE(eacP1) || MODE(eacP2) || MODE(eacP3))
 181             {
 182                 unsigned int mmm;
 183
 184                 for (m = 0; (m < DIM); m++)
 185                 {
 186                     xj[m] = c1[j3 + m];
 187                     xk[m] = c1[jk3 + m];
 188                 }
 189                 cth = cos_angle(xj, xk);
 190
 191                 if (cth - 1.0 > 1.0e-15)
 192                 {
 193                     printf("j: %d, k: %d, xj:(%g,%g,%g), xk:(%g,%g,%g)\n",
 194                            j,
 195                            k,
 196                            xj[XX],
 197                            xj[YY],
 198                            xj[ZZ],
 199                            xk[XX],
 200                            xk[YY],
 201                            xk[ZZ]);
 202                 }
 203                 mmm = 1;
 204                 if (MODE(eacP2))
 205                 {
 206                     mmm = 2;
 207                 }
 208                 else if (MODE(eacP3))
 209                 {
 210                     mmm = 3;
 211                 }
 212                 corr[k] += LegendreP(cth, mmm); /* 1.5*cth*cth-0.5; */
 213             }
 214             else if (MODE(eacRcross))
 215             {
 216                 rvec xj, xk, rr;
 217                 for (m = 0; (m < DIM); m++)
 218                 {
 219                     xj[m] = c1[j3 + m];
 220                     xk[m] = c1[jk3 + m];
 221                 }
 222                 cprod(xj, xk, rr);
 223
 224                 corr[k] += iprod(rr, rr);
 225             }
 226             else if (MODE(eacVector))
 227             {
 228                 for (m = 0; (m < DIM); m++)
 229                 {
 230                     xj[m] = c1[j3 + m];
 231                     xk[m] = c1[jk3 + m];
 232                 }
 233                 ccc = iprod(xj, xk);
 234
 235                 corr[k] += ccc;
 236             }
 237             else
 238             {
 239                 gmx_fatal(FARGS, "\nInvalid mode (%lu) in do_ac_core", mode);
 240             }
 241         }
 242     }
 243     /* Correct for the number of points and copy results to the data array */
 244     for (j = 0; (j < nout); j++)
 245     {
 246         n     = (nframes - j + (nrestart - 1)) / nrestart;
 247         c1[j] = corr[j] / n;
 248     }
 249 }
 250
 251 /*! \brief Routine to normalize ACF, dividing by corr[0]. */
 252 static void normalize_acf(int nout, real corr[])
 253 {
 254     int    j;
 255     double c0;
 256
 257     if (debug)
 258     {
 259         fprintf(debug, "Before normalization\n");
 260         for (j = 0; (j < nout); j++)
 261         {
 262             fprintf(debug, "%5d  %10f\n", j, corr[j]);
 263         }
 264     }
 265
 266     /* Normalisation makes that c[0] = 1.0 and that other points are scaled
 267      * accordingly.
 268      */
 269     if (fabs(corr[0]) < 1e-5)
 270     {
 271         c0 = 1.0;
 272     }
 273     else
 274     {
 275         c0 = 1.0 / corr[0];
 276     }
 277     for (j = 0; (j < nout); j++)
 278     {
 279         corr[j] *= c0;
 280     }
 281
 282     if (debug)
 283     {
 284         fprintf(debug, "After normalization\n");
 285         for (j = 0; (j < nout); j++)
 286         {
 287             fprintf(debug, "%5d  %10f\n", j, corr[j]);
 288         }
 289     }
 290 }
 291
 292 /*! \brief Routine that averages ACFs. */
 293 static void average_acf(gmx_bool bVerbose, int n, int nitem, real** c1)
 294 {
 295     real c0;
 296     int  i, j;
 297
 298     if (bVerbose)
 299     {
 300         printf("Averaging correlation functions\n");
 301     }
 302
 303     for (j = 0; (j < n); j++)
 304     {
 305         c0 = 0;
 306         for (i = 0; (i < nitem); i++)
 307         {
 308             c0 += c1[i][j];
 309         }
 310         c1[0][j] = c0 / nitem;
 311     }
 312 }
 313
 314 /*! \brief Normalize ACFs. */
 315 static void norm_and_scale_vectors(int nframes, real c1[], real scale)
 316 {
 317     int   j, m;
 318     real* rij;
 319
 320     for (j = 0; (j < nframes); j++)
 321     {
 322         rij = &(c1[j * DIM]);
 323         unitv(rij, rij);
 324         for (m = 0; (m < DIM); m++)
 325         {
 326             rij[m] *= scale;
 327         }
 328     }
 329 }
 330
 331 /*! \brief Debugging */
 332 static void dump_tmp(char* s, int n, real c[])
 333 {
 334     FILE* fp;
 335     int   i;
 336
 337     fp = gmx_ffopen(s, "w");
 338     for (i = 0; (i < n); i++)
 339     {
 340         fprintf(fp, "%10d  %10g\n", i, c[i]);
 341     }
 342     gmx_ffclose(fp);
 343 }
 344
 345 /*! \brief High level ACF routine. */
 346 static void do_four_core(unsigned long mode, int nframes, real c1[], real csum[], real ctmp[])
 347 {
 348     real* cfour;
 349     char  buf[32];
 350     real  fac;
 351     int   j, m, m1;
 352
 353     snew(cfour, nframes);
 354
 355     if (MODE(eacNormal))
 356     {
 357         /********************************************
 358          *  N O R M A L
 359          ********************************************/
 360         low_do_four_core(nframes, c1, csum, enNorm);
 361     }
 362     else if (MODE(eacCos))
 363     {
 364         /***************************************************
 365          * C O S I N E
 366          ***************************************************/
 367         /* Copy the data to temp array. Since we need it twice
 368          * we can't overwrite original.
 369          */
 370         for (j = 0; (j < nframes); j++)
 371         {
 372             ctmp[j] = c1[j];
 373         }
 374
 375         /* Cosine term of AC function */
 376         low_do_four_core(nframes, ctmp, cfour, enCos);
 377         for (j = 0; (j < nframes); j++)
 378         {
 379             c1[j] = cfour[j];
 380         }
 381
 382         /* Sine term of AC function */
 383         low_do_four_core(nframes, ctmp, cfour, enSin);
 384         for (j = 0; (j < nframes); j++)
 385         {
 386             c1[j] += cfour[j];
 387             csum[j] = c1[j];
 388         }
 389     }
 390     else if (MODE(eacP2))
 391     {
 392         /***************************************************
 393          * Legendre polynomials
 394          ***************************************************/
 395         /* First normalize the vectors */
 396         norm_and_scale_vectors(nframes, c1, 1.0);
 397
 398         /* For P2 thingies we have to do six FFT based correls
 399          * First for XX^2, then for YY^2, then for ZZ^2
 400          * Then we have to do XY, YZ and XZ (counting these twice)
 401          * After that we sum them and normalise
 402          * P2(x) = (3 * cos^2 (x) - 1)/2
 403          * for unit vectors u and v we compute the cosine as the inner product
 404          * cos(u,v) = uX vX + uY vY + uZ vZ
 405          *
 406          *        oo
 407          *        /
 408          * C(t) = |  (3 cos^2(u(t'),u(t'+t)) - 1)/2 dt'
 409          *        /
 410          *        0
 411          *
 412          * For ACF we need:
 413          * P2(u(0),u(t)) = [3 * (uX(0) uX(t) +
 414          *                       uY(0) uY(t) +
 415          *                       uZ(0) uZ(t))^2 - 1]/2
 416          *               = [3 * ((uX(0) uX(t))^2 +
 417          *                       (uY(0) uY(t))^2 +
 418          *                       (uZ(0) uZ(t))^2 +
 419          *                 2(uX(0) uY(0) uX(t) uY(t)) +
 420          *                 2(uX(0) uZ(0) uX(t) uZ(t)) +
 421          *                 2(uY(0) uZ(0) uY(t) uZ(t))) - 1]/2
 422          *
 423          *               = [(3/2) * (<uX^2> + <uY^2> + <uZ^2> +
 424          *                         2<uXuY> + 2<uXuZ> + 2<uYuZ>) - 0.5]
 425          *
 426          */
 427
 428         /* Because of normalization the number of -0.5 to subtract
 429          * depends on the number of data points!
 430          */
 431         for (j = 0; (j < nframes); j++)
 432         {
 433             csum[j] = -0.5 * (nframes - j);
 434         }
 435
 436         /***** DIAGONAL ELEMENTS ************/
 437         for (m = 0; (m < DIM); m++)
 438         {
 439             /* Copy the vector data in a linear array */
 440             for (j = 0; (j < nframes); j++)
 441             {
 442                 ctmp[j] = gmx::square(c1[DIM * j + m]);
 443             }
 444             if (debug)
 445             {
 446                 sprintf(buf, "c1diag%d.xvg", m);
 447                 dump_tmp(buf, nframes, ctmp);
 448             }
 449
 450             low_do_four_core(nframes, ctmp, cfour, enNorm);
 451
 452             if (debug)
 453             {
 454                 sprintf(buf, "c1dfout%d.xvg", m);
 455                 dump_tmp(buf, nframes, cfour);
 456             }
 457             fac = 1.5;
 458             for (j = 0; (j < nframes); j++)
 459             {
 460                 csum[j] += fac * (cfour[j]);
 461             }
 462         }
 463         /******* OFF-DIAGONAL ELEMENTS **********/
 464         for (m = 0; (m < DIM); m++)
 465         {
 466             /* Copy the vector data in a linear array */
 467             m1 = (m + 1) % DIM;
 468             for (j = 0; (j < nframes); j++)
 469             {
 470                 ctmp[j] = c1[DIM * j + m] * c1[DIM * j + m1];
 471             }
 472
 473             if (debug)
 474             {
 475                 sprintf(buf, "c1off%d.xvg", m);
 476                 dump_tmp(buf, nframes, ctmp);
 477             }
 478             low_do_four_core(nframes, ctmp, cfour, enNorm);
 479             if (debug)
 480             {
 481                 sprintf(buf, "c1ofout%d.xvg", m);
 482                 dump_tmp(buf, nframes, cfour);
 483             }
 484             fac = 3.0;
 485             for (j = 0; (j < nframes); j++)
 486             {
 487                 csum[j] += fac * cfour[j];
 488             }
 489         }
 490     }
 491     else if (MODE(eacP1) || MODE(eacVector))
 492     {
 493         /***************************************************
 494          * V E C T O R & P1
 495          ***************************************************/
 496         if (MODE(eacP1))
 497         {
 498             /* First normalize the vectors */
 499             norm_and_scale_vectors(nframes, c1, 1.0);
 500         }
 501
 502         /* For vector thingies we have to do three FFT based correls
 503          * First for XX, then for YY, then for ZZ
 504          * After that we sum them and normalise
 505          */
 506         for (j = 0; (j < nframes); j++)
 507         {
 508             csum[j] = 0.0;
 509         }
 510         for (m = 0; (m < DIM); m++)
 511         {
 512             /* Copy the vector data in a linear array */
 513             for (j = 0; (j < nframes); j++)
 514             {
 515                 ctmp[j] = c1[DIM * j + m];
 516             }
 517             low_do_four_core(nframes, ctmp, cfour, enNorm);
 518             for (j = 0; (j < nframes); j++)
 519             {
 520                 csum[j] += cfour[j];
 521             }
 522         }
 523     }
 524     else
 525     {
 526         gmx_fatal(FARGS, "\nUnknown mode in do_autocorr (%lu)", mode);
 527     }
 528
 529     sfree(cfour);
 530     for (j = 0; (j < nframes); j++)
 531     {
 532         c1[j] = csum[j] / static_cast<real>(nframes - j);
 533     }
 534 }
 535
 536 void low_do_autocorr(const char*             fn,
 537                      const gmx_output_env_t* oenv,
 538                      const char*             title,
 539                      int                     nframes,
 540                      int                     nitem,
 541                      int                     nout,
 542                      real**                  c1,
 543                      real                    dt,
 544                      unsigned long           mode,
 545                      int                     nrestart,
 546                      gmx_bool                bAver,
 547                      gmx_bool                bNormalize,
 548                      gmx_bool                bVerbose,
 549                      real                    tbeginfit,
 550                      real                    tendfit,
 551                      int                     eFitFn)
 552 {
 553     FILE *   fp, *gp = nullptr;
 554     int      i;
 555     real*    csum;
 556     real *   ctmp, *fit;
 557     real     sum, Ct2av, Ctav;
 558     gmx_bool bFour = acf.bFour;
 559
 560     /* Check flags and parameters */
 561     nout = get_acfnout();
 562     if (nout == -1)
 563     {
 564         nout = acf.nout = (nframes + 1) / 2;
 565     }
 566     else if (nout > nframes)
 567     {
 568         nout = nframes;
 569     }
 570
 571     if (MODE(eacCos) && MODE(eacVector))
 572     {
 573         gmx_fatal(FARGS, "Incompatible options bCos && bVector (%s, %d)", __FILE__, __LINE__);
 574     }
 575     if ((MODE(eacP3) || MODE(eacRcross)) && bFour)
 576     {
 577         if (bVerbose)
 578         {
 579             fprintf(stderr, "Can't combine mode %lu with FFT, turning off FFT\n", mode);
 580         }
 581         bFour = FALSE;
 582     }
 583     if (MODE(eacNormal) && MODE(eacVector))
 584     {
 585         gmx_fatal(FARGS, "Incompatible mode bits: normal and vector (or Legendre)");
 586     }
 587
 588     /* Print flags and parameters */
 589     if (bVerbose)
 590     {
 591         printf("Will calculate %s of %d thingies for %d frames\n", title ? title : "autocorrelation", nitem, nframes);
 592         printf("bAver = %s, bFour = %s bNormalize= %s\n",
 593                gmx::boolToString(bAver),
 594                gmx::boolToString(bFour),
 595                gmx::boolToString(bNormalize));
 596         printf("mode = %lu, dt = %g, nrestart = %d\n", mode, dt, nrestart);
 597     }
 598     /* Allocate temp arrays */
 599     snew(csum, nframes);
 600     snew(ctmp, nframes);
 601
 602     /* Loop over items (e.g. molecules or dihedrals)
 603      * In this loop the actual correlation functions are computed, but without
 604      * normalizing them.
 605      */
 606     for (int i = 0; i < nitem; i++)
 607     {
 608         if (bVerbose && (((i % 100) == 0) || (i == nitem - 1)))
 609         {
 610             fprintf(stderr, "\rThingie %d", i + 1);
 611             fflush(stderr);
 612         }
 613
 614         if (bFour)
 615         {
 616             do_four_core(mode, nframes, c1[i], csum, ctmp);
 617         }
 618         else
 619         {
 620             do_ac_core(nframes, nout, ctmp, c1[i], nrestart, mode);
 621         }
 622     }
 623     if (bVerbose)
 624     {
 625         fprintf(stderr, "\n");
 626     }
 627     sfree(ctmp);
 628     sfree(csum);
 629
 630     if (fn)
 631     {
 632         snew(fit, nout);
 633         fp = xvgropen(fn, title, "Time (ps)", "C(t)", oenv);
 634     }
 635     else
 636     {
 637         fit = nullptr;
 638         fp  = nullptr;
 639     }
 640     if (bAver)
 641     {
 642         if (nitem > 1)
 643         {
 644             average_acf(bVerbose, nframes, nitem, c1);
 645         }
 646
 647         if (bNormalize)
 648         {
 649             normalize_acf(nout, c1[0]);
 650         }
 651
 652         if (eFitFn != effnNONE)
 653         {
 654             fit_acf(nout, eFitFn, oenv, fn != nullptr, tbeginfit, tendfit, dt, c1[0], fit);
 655             sum = print_and_integrate(fp, nout, dt, c1[0], fit, 1);
 656         }
 657         else
 658         {
 659             sum = print_and_integrate(fp, nout, dt, c1[0], nullptr, 1);
 660         }
 661         if (bVerbose)
 662         {
 663             printf("Correlation time (integral over corrfn): %g (ps)\n", sum);
 664         }
 665     }
 666     else
 667     {
 668         /* Not averaging. Normalize individual ACFs */
 669         Ctav = Ct2av = 0;
 670         if (debug)
 671         {
 672             gp = xvgropen("ct-distr.xvg", "Correlation times", "item", "time (ps)", oenv);
 673         }
 674         for (i = 0; i < nitem; i++)
 675         {
 676             if (bNormalize)
 677             {
 678                 normalize_acf(nout, c1[i]);
 679             }
 680             if (eFitFn != effnNONE)
 681             {
 682                 fit_acf(nout, eFitFn, oenv, fn != nullptr, tbeginfit, tendfit, dt, c1[i], fit);
 683                 sum = print_and_integrate(fp, nout, dt, c1[i], fit, 1);
 684             }
 685             else
 686             {
 687                 sum = print_and_integrate(fp, nout, dt, c1[i], nullptr, 1);
 688                 if (debug)
 689                 {
 690                     fprintf(debug, "CORRelation time (integral over corrfn %d): %g (ps)\n", i, sum);
 691                 }
 692             }
 693             Ctav += sum;
 694             Ct2av += sum * sum;
 695             if (debug)
 696             {
 697                 fprintf(gp, "%5d  %.3f\n", i, sum);
 698             }
 699         }
 700         if (debug)
 701         {
 702             xvgrclose(gp);
 703         }
 704         if (nitem > 1)
 705         {
 706             Ctav /= nitem;
 707             Ct2av /= nitem;
 708             printf("Average correlation time %.3f Std. Dev. %.3f Error %.3f (ps)\n",
 709                    Ctav,
 710                    std::sqrt((Ct2av - gmx::square(Ctav))),
 711                    std::sqrt((Ct2av - gmx::square(Ctav)) / (nitem - 1)));
 712         }
 713     }
 714     if (fp)
 715     {
 716         xvgrclose(fp);
 717     }
 718     sfree(fit);
 719 }
 720
 721 /*! \brief Legend for selecting Legendre polynomials. */
 722 static const char* Leg[] = { nullptr, "0", "1", "2", "3", nullptr };
 723
 724 t_pargs* add_acf_pargs(int* npargs, t_pargs* pa)
 725 {
 726     t_pargs acfpa[] = {
 727         { "-acflen",
 728           FALSE,
 729           etINT,
 730           { &acf.nout },
 731           "Length of the ACF, default is half the number of frames" },
 732         { "-normalize", FALSE, etBOOL, { &acf.bNormalize }, "Normalize ACF" },
 733         { "-fftcorr",
 734           FALSE,
 735           etBOOL,
 736           { &acf.bFour },
 737           "HIDDENUse fast fourier transform for correlation function" },
 738         { "-nrestart",
 739           FALSE,
 740           etINT,
 741           { &acf.nrestart },
 742           "HIDDENNumber of frames between time origins for ACF when no FFT is used" },
 743         { "-P", FALSE, etENUM, { Leg }, "Order of Legendre polynomial for ACF (0 indicates none)" },
 744         { "-fitfn", FALSE, etENUM, { s_ffn }, "Fit function" },
 745         { "-beginfit",
 746           FALSE,
 747           etREAL,
 748           { &acf.tbeginfit },
 749           "Time where to begin the exponential fit of the correlation function" },
 750         { "-endfit",
 751           FALSE,
 752           etREAL,
 753           { &acf.tendfit },
 754           "Time where to end the exponential fit of the correlation function, -1 is until the "
 755           "end" },
 756     };
 757     t_pargs* ppa;
 758     int      i, npa;
 759
 760     npa = asize(acfpa);
 761     snew(ppa, *npargs + npa);
 762     for (i = 0; (i < *npargs); i++)
 763     {
 764         ppa[i] = pa[i];
 765     }
 766     for (i = 0; (i < npa); i++)
 767     {
 768         ppa[*npargs + i] = acfpa[i];
 769     }
 770     (*npargs) += npa;
 771
 772     acf.mode       = 0;
 773     acf.nrestart   = 1;
 774     acf.nout       = -1;
 775     acf.P          = 0;
 776     acf.fitfn      = effnEXP1;
 777     acf.bFour      = TRUE;
 778     acf.bNormalize = TRUE;
 779     acf.tbeginfit  = 0.0;
 780     acf.tendfit    = -1;
 781
 782     bACFinit = TRUE;
 783
 784     return ppa;
 785 }
 786
 787 void do_autocorr(const char*             fn,
 788                  const gmx_output_env_t* oenv,
 789                  const char*             title,
 790                  int                     nframes,
 791                  int                     nitem,
 792                  real**                  c1,
 793                  real                    dt,
 794                  unsigned long           mode,
 795                  gmx_bool                bAver)
 796 {
 797     if (!bACFinit)
 798     {
 799         printf("ACF data structures have not been initialised. Call add_acf_pargs\n");
 800     }
 801
 802     /* Handle enumerated types */
 803     sscanf(Leg[0], "%d", &acf.P);
 804     acf.fitfn = sffn2effn(s_ffn);
 805
 806     switch (acf.P)
 807     {
 808         case 1: mode = mode | eacP1; break;
 809         case 2: mode = mode | eacP2; break;
 810         case 3: mode = mode | eacP3; break;
 811         default: break;
 812     }
 813
 814     low_do_autocorr(fn,
 815                     oenv,
 816                     title,
 817                     nframes,
 818                     nitem,
 819                     acf.nout,
 820                     c1,
 821                     dt,
 822                     mode,
 823                     acf.nrestart,
 824                     bAver,
 825                     acf.bNormalize,
 826                     bDebugMode(),
 827                     acf.tbeginfit,
 828                     acf.tendfit,
 829                     acf.fitfn);
 830 }
 831
 832 int get_acfnout()
 833 {
 834     if (!bACFinit)
 835     {
 836         gmx_fatal(FARGS, "ACF data not initialized yet");
 837     }
 838
 839     return acf.nout;
 840 }
 841
 842 int get_acffitfn()
 843 {
 844     if (!bACFinit)
 845     {
 846         gmx_fatal(FARGS, "ACF data not initialized yet");
 847     }
 848
 849     return sffn2effn(s_ffn);
 850 }