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