src/gromacs/correlationfunctions/autocorr.c

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