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