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