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