src/gromacs/gmxana/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, 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 #include "gmxpre.h"
  38
  39 #include "config.h"
  40
  41 #include <math.h>
  42 #include <stdio.h>
  43 #include <string.h>
  44
  45 #include "gromacs/legacyheaders/macros.h"
  46 #include "gromacs/legacyheaders/typedefs.h"
  47 #include "gromacs/utility/smalloc.h"
  48 #include "gromacs/fileio/xvgr.h"
  49 #include "gromacs/utility/futil.h"
  50 #include "gstat.h"
  51 #include "gromacs/legacyheaders/names.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/math/vec.h"
  54 #include "correl.h"
  55
  56 #define MODE(x) ((mode & (x)) == (x))
  57
  58 typedef struct {
  59     unsigned long mode;
  60     int           nrestart, nout, P, fitfn;
  61     gmx_bool      bFour, bNormalize;
  62     real          tbeginfit, tendfit;
  63 } t_acf;
  64
  65 static gmx_bool  bACFinit = FALSE;
  66 static t_acf     acf;
  67
  68 enum {
  69     enNorm, enCos, enSin
  70 };
  71
  72 int sffn2effn(const char **sffn)
  73 {
  74     int eFitFn, i;
  75
  76     eFitFn = 0;
  77     for (i = 0; i < effnNR; i++)
  78     {
  79         if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0)
  80         {
  81             eFitFn = i;
  82         }
  83     }
  84
  85     return eFitFn;
  86 }
  87
  88 static void low_do_four_core(int nfour, int nframes, real c1[], real cfour[],
  89                              int nCos, gmx_bool bPadding)
  90 {
  91     int  i = 0;
  92     real aver, *ans;
  93
  94     aver = 0.0;
  95     switch (nCos)
  96     {
  97         case enNorm:
  98             for (i = 0; (i < nframes); i++)
  99             {
 100                 aver    += c1[i];
 101                 cfour[i] = c1[i];
 102             }
 103             break;
 104         case enCos:
 105             for (i = 0; (i < nframes); i++)
 106             {
 107                 cfour[i] = cos(c1[i]);
 108             }
 109             break;
 110         case enSin:
 111             for (i = 0; (i < nframes); i++)
 112             {
 113                 cfour[i] = sin(c1[i]);
 114             }
 115             break;
 116         default:
 117             gmx_fatal(FARGS, "nCos = %d, %s %d", nCos, __FILE__, __LINE__);
 118     }
 119     for (; (i < nfour); i++)
 120     {
 121         cfour[i] = 0.0;
 122     }
 123
 124     if (bPadding)
 125     {
 126         aver /= nframes;
 127         /* printf("aver = %g\n",aver); */
 128         for (i = 0; (i < nframes); i++)
 129         {
 130             cfour[i] -= aver;
 131         }
 132     }
 133
 134     snew(ans, 2*nfour);
 135     correl(cfour-1, cfour-1, nfour, ans-1);
 136
 137     if (bPadding)
 138     {
 139         for (i = 0; (i < nfour); i++)
 140         {
 141             cfour[i] = ans[i]+sqr(aver);
 142         }
 143     }
 144     else
 145     {
 146         for (i = 0; (i < nfour); i++)
 147         {
 148             cfour[i] = ans[i];
 149         }
 150     }
 151
 152     sfree(ans);
 153 }
 154
 155 static void do_ac_core(int nframes, int nout,
 156                        real corr[], real c1[], int nrestart,
 157                        unsigned long mode)
 158 {
 159     int     j, k, j3, jk3, m, n;
 160     real    ccc, cth;
 161     rvec    xj, xk, rr;
 162
 163     if (nrestart < 1)
 164     {
 165         printf("WARNING: setting number of restarts to 1\n");
 166         nrestart = 1;
 167     }
 168     if (debug)
 169     {
 170         fprintf(debug,
 171                 "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n",
 172                 nframes, nout, nrestart, mode);
 173     }
 174
 175     for (j = 0; (j < nout); j++)
 176     {
 177         corr[j] = 0;
 178     }
 179
 180     /* Loop over starting points. */
 181     for (j = 0; (j < nframes); j += nrestart)
 182     {
 183         j3  = DIM*j;
 184
 185         /* Loop over the correlation length for this starting point */
 186         for (k = 0; (k < nout) && (j+k < nframes); k++)
 187         {
 188             jk3 = DIM*(j+k);
 189
 190             /* Switch over possible ACF types.
 191              * It might be more efficient to put the loops inside the switch,
 192              * but this is more clear, and save development time!
 193              */
 194             if (MODE(eacNormal))
 195             {
 196                 corr[k] += c1[j]*c1[j+k];
 197             }
 198             else if (MODE(eacCos))
 199             {
 200                 /* Compute the cos (phi(t)-phi(t+dt)) */
 201                 corr[k] += cos(c1[j]-c1[j+k]);
 202             }
 203             else if (MODE(eacIden))
 204             {
 205                 /* Check equality (phi(t)==phi(t+dt)) */
 206                 corr[k] += (c1[j] == c1[j+k]) ? 1 : 0;
 207             }
 208             else if (MODE(eacP1) || MODE(eacP2) || MODE(eacP3))
 209             {
 210                 for (m = 0; (m < DIM); m++)
 211                 {
 212                     xj[m] = c1[j3+m];
 213                     xk[m] = c1[jk3+m];
 214                 }
 215                 cth = cos_angle(xj, xk);
 216
 217                 if (cth-1.0 > 1.0e-15)
 218                 {
 219                     printf("j: %d, k: %d, xj:(%g,%g,%g), xk:(%g,%g,%g)\n",
 220                            j, k, xj[XX], xj[YY], xj[ZZ], xk[XX], xk[YY], xk[ZZ]);
 221                 }
 222
 223                 corr[k] += LegendreP(cth, mode); /* 1.5*cth*cth-0.5; */
 224             }
 225             else if (MODE(eacRcross))
 226             {
 227                 for (m = 0; (m < DIM); m++)
 228                 {
 229                     xj[m] = c1[j3+m];
 230                     xk[m] = c1[jk3+m];
 231                 }
 232                 cprod(xj, xk, rr);
 233
 234                 corr[k] += iprod(rr, rr);
 235             }
 236             else if (MODE(eacVector))
 237             {
 238                 for (m = 0; (m < DIM); m++)
 239                 {
 240                     xj[m] = c1[j3+m];
 241                     xk[m] = c1[jk3+m];
 242                 }
 243                 ccc = iprod(xj, xk);
 244
 245                 corr[k] += ccc;
 246             }
 247             else
 248             {
 249                 gmx_fatal(FARGS, "\nInvalid mode (%d) in do_ac_core", mode);
 250             }
 251         }
 252     }
 253     /* Correct for the number of points and copy results to the data array */
 254     for (j = 0; (j < nout); j++)
 255     {
 256         n     = (nframes-j+(nrestart-1))/nrestart;
 257         c1[j] = corr[j]/n;
 258     }
 259 }
 260
 261 void normalize_acf(int nout, real corr[])
 262 {
 263     int  j;
 264     real c0;
 265
 266     if (debug)
 267     {
 268         fprintf(debug, "Before normalization\n");
 269         for (j = 0; (j < nout); j++)
 270         {
 271             fprintf(debug, "%5d  %10f\n", j, corr[j]);
 272         }
 273     }
 274
 275     /* Normalisation makes that c[0] = 1.0 and that other points are scaled
 276      * accordingly.
 277      */
 278     if (corr[0] == 0.0)
 279     {
 280         c0 = 1.0;
 281     }
 282     else
 283     {
 284         c0 = 1.0/corr[0];
 285     }
 286     for (j = 0; (j < nout); j++)
 287     {
 288         corr[j] *= c0;
 289     }
 290
 291     if (debug)
 292     {
 293         fprintf(debug, "After normalization\n");
 294         for (j = 0; (j < nout); j++)
 295         {
 296             fprintf(debug, "%5d  %10f\n", j, corr[j]);
 297         }
 298     }
 299 }
 300
 301 void average_acf(gmx_bool bVerbose, int n, int nitem, real **c1)
 302 {
 303     real c0;
 304     int  i, j;
 305
 306     if (bVerbose)
 307     {
 308         printf("Averaging correlation functions\n");
 309     }
 310
 311     for (j = 0; (j < n); j++)
 312     {
 313         c0 = 0;
 314         for (i = 0; (i < nitem); i++)
 315         {
 316             c0 += c1[i][j];
 317         }
 318         c1[0][j] = c0/nitem;
 319     }
 320 }
 321
 322 void norm_and_scale_vectors(int nframes, real c1[], real scale)
 323 {
 324     int   j, m;
 325     real *rij;
 326
 327     for (j = 0; (j < nframes); j++)
 328     {
 329         rij = &(c1[j*DIM]);
 330         unitv(rij, rij);
 331         for (m = 0; (m < DIM); m++)
 332         {
 333             rij[m] *= scale;
 334         }
 335     }
 336 }
 337
 338 void dump_tmp(char *s, int n, real c[])
 339 {
 340     FILE *fp;
 341     int   i;
 342
 343     fp = gmx_ffopen(s, "w");
 344     for (i = 0; (i < n); i++)
 345     {
 346         fprintf(fp, "%10d  %10g\n", i, c[i]);
 347     }
 348     gmx_ffclose(fp);
 349 }
 350
 351 real print_and_integrate(FILE *fp, int n, real dt, real c[], real *fit, int nskip)
 352 {
 353     real c0, sum;
 354     int  j;
 355
 356     /* Use trapezoidal rule for calculating integral */
 357     sum = 0.0;
 358     for (j = 0; (j < n); j++)
 359     {
 360         c0 = c[j];
 361         if (fp && (nskip == 0 || j % nskip == 0))
 362         {
 363             fprintf(fp, "%10.3f  %10.5f\n", j*dt, c0);
 364         }
 365         if (j > 0)
 366         {
 367             sum += dt*(c0+c[j-1]);
 368         }
 369     }
 370     if (fp)
 371     {
 372         fprintf(fp, "&\n");
 373         if (fit)
 374         {
 375             for (j = 0; (j < n); j++)
 376             {
 377                 if (nskip == 0 || j % nskip == 0)
 378                 {
 379                     fprintf(fp, "%10.3f  %10.5f\n", j*dt, fit[j]);
 380                 }
 381             }
 382             fprintf(fp, "&\n");
 383         }
 384     }
 385     return sum*0.5;
 386 }
 387
 388 real evaluate_integral(int n, real x[], real y[], real dy[], real aver_start,
 389                        real *stddev)
 390 {
 391     double sum, sum_var, w;
 392     double sum_tail = 0, sum2_tail = 0;
 393     int    j, nsum_tail = 0;
 394
 395     /* Use trapezoidal rule for calculating integral */
 396     if (n <= 0)
 397     {
 398         gmx_fatal(FARGS, "Evaluating integral: n = %d (file %s, line %d)",
 399                   n, __FILE__, __LINE__);
 400     }
 401
 402     sum     = 0;
 403     sum_var = 0;
 404     for (j = 0; (j < n); j++)
 405     {
 406         w = 0;
 407         if (j > 0)
 408         {
 409             w += 0.5*(x[j] - x[j-1]);
 410         }
 411         if (j < n-1)
 412         {
 413             w += 0.5*(x[j+1] - x[j]);
 414         }
 415         sum += w*y[j];
 416         if (dy)
 417         {
 418             /* Assume all errors are uncorrelated */
 419             sum_var += sqr(w*dy[j]);
 420         }
 421
 422         if ((aver_start > 0) && (x[j] >= aver_start))
 423         {
 424             sum_tail  += sum;
 425             sum2_tail += sqrt(sum_var);
 426             nsum_tail += 1;
 427         }
 428     }
 429
 430     if (nsum_tail > 0)
 431     {
 432         sum = sum_tail/nsum_tail;
 433         /* This is a worst case estimate, assuming all stddev's are correlated. */
 434         *stddev = sum2_tail/nsum_tail;
 435     }
 436     else
 437     {
 438         *stddev = sqrt(sum_var);
 439     }
 440
 441     return sum;
 442 }
 443
 444 void do_four_core(unsigned long mode, int nfour, int nf2, int nframes,
 445                   real c1[], real csum[], real ctmp[])
 446 {
 447     real   *cfour;
 448     char    buf[32];
 449     real    fac;
 450     int     j, m, m1;
 451
 452     snew(cfour, nfour);
 453
 454     if (MODE(eacNormal))
 455     {
 456         /********************************************
 457          *  N O R M A L
 458          ********************************************/
 459         low_do_four_core(nfour, nf2, c1, csum, enNorm, FALSE);
 460     }
 461     else if (MODE(eacCos))
 462     {
 463         /***************************************************
 464          * C O S I N E
 465          ***************************************************/
 466         /* Copy the data to temp array. Since we need it twice
 467          * we can't overwrite original.
 468          */
 469         for (j = 0; (j < nf2); j++)
 470         {
 471             ctmp[j] = c1[j];
 472         }
 473
 474         /* Cosine term of AC function */
 475         low_do_four_core(nfour, nf2, ctmp, cfour, enCos, FALSE);
 476         for (j = 0; (j < nf2); j++)
 477         {
 478             c1[j]  = cfour[j];
 479         }
 480
 481         /* Sine term of AC function */
 482         low_do_four_core(nfour, nf2, ctmp, cfour, enSin, FALSE);
 483         for (j = 0; (j < nf2); j++)
 484         {
 485             c1[j]  += cfour[j];
 486             csum[j] = c1[j];
 487         }
 488     }
 489     else if (MODE(eacP2))
 490     {
 491         /***************************************************
 492          * Legendre polynomials
 493          ***************************************************/
 494         /* First normalize the vectors */
 495         norm_and_scale_vectors(nframes, c1, 1.0);
 496
 497         /* For P2 thingies we have to do six FFT based correls
 498          * First for XX^2, then for YY^2, then for ZZ^2
 499          * Then we have to do XY, YZ and XZ (counting these twice)
 500          * After that we sum them and normalise
 501          * P2(x) = (3 * cos^2 (x) - 1)/2
 502          * for unit vectors u and v we compute the cosine as the inner product
 503          * cos(u,v) = uX vX + uY vY + uZ vZ
 504          *
 505          *        oo
 506          *        /
 507          * C(t) = |  (3 cos^2(u(t'),u(t'+t)) - 1)/2 dt'
 508          *        /
 509          *        0
 510          *
 511          * For ACF we need:
 512          * P2(u(0),u(t)) = [3 * (uX(0) uX(t) +
 513          *                       uY(0) uY(t) +
 514          *                       uZ(0) uZ(t))^2 - 1]/2
 515          *               = [3 * ((uX(0) uX(t))^2 +
 516          *                       (uY(0) uY(t))^2 +
 517          *                       (uZ(0) uZ(t))^2 +
 518          *                 2(uX(0) uY(0) uX(t) uY(t)) +
 519          *                 2(uX(0) uZ(0) uX(t) uZ(t)) +
 520          *                 2(uY(0) uZ(0) uY(t) uZ(t))) - 1]/2
 521          *
 522          *               = [(3/2) * (<uX^2> + <uY^2> + <uZ^2> +
 523          *                         2<uXuY> + 2<uXuZ> + 2<uYuZ>) - 0.5]
 524          *
 525          */
 526
 527         /* Because of normalization the number of -0.5 to subtract
 528          * depends on the number of data points!
 529          */
 530         for (j = 0; (j < nf2); j++)
 531         {
 532             csum[j]  = -0.5*(nf2-j);
 533         }
 534
 535         /***** DIAGONAL ELEMENTS ************/
 536         for (m = 0; (m < DIM); m++)
 537         {
 538             /* Copy the vector data in a linear array */
 539             for (j = 0; (j < nf2); j++)
 540             {
 541                 ctmp[j]  = sqr(c1[DIM*j+m]);
 542             }
 543             if (debug)
 544             {
 545                 sprintf(buf, "c1diag%d.xvg", m);
 546                 dump_tmp(buf, nf2, ctmp);
 547             }
 548
 549             low_do_four_core(nfour, nf2, ctmp, cfour, enNorm, FALSE);
 550
 551             if (debug)
 552             {
 553                 sprintf(buf, "c1dfout%d.xvg", m);
 554                 dump_tmp(buf, nf2, cfour);
 555             }
 556             fac = 1.5;
 557             for (j = 0; (j < nf2); j++)
 558             {
 559                 csum[j] += fac*(cfour[j]);
 560             }
 561         }
 562         /******* OFF-DIAGONAL ELEMENTS **********/
 563         for (m = 0; (m < DIM); m++)
 564         {
 565             /* Copy the vector data in a linear array */
 566             m1 = (m+1) % DIM;
 567             for (j = 0; (j < nf2); j++)
 568             {
 569                 ctmp[j] = c1[DIM*j+m]*c1[DIM*j+m1];
 570             }
 571
 572             if (debug)
 573             {
 574                 sprintf(buf, "c1off%d.xvg", m);
 575                 dump_tmp(buf, nf2, ctmp);
 576             }
 577             low_do_four_core(nfour, nf2, ctmp, cfour, enNorm, FALSE);
 578             if (debug)
 579             {
 580                 sprintf(buf, "c1ofout%d.xvg", m);
 581                 dump_tmp(buf, nf2, cfour);
 582             }
 583             fac = 3.0;
 584             for (j = 0; (j < nf2); j++)
 585             {
 586                 csum[j] += fac*cfour[j];
 587             }
 588         }
 589     }
 590     else if (MODE(eacP1) || MODE(eacVector))
 591     {
 592         /***************************************************
 593          * V E C T O R & P1
 594          ***************************************************/
 595         if (MODE(eacP1))
 596         {
 597             /* First normalize the vectors */
 598             norm_and_scale_vectors(nframes, c1, 1.0);
 599         }
 600
 601         /* For vector thingies we have to do three FFT based correls
 602          * First for XX, then for YY, then for ZZ
 603          * After that we sum them and normalise
 604          */
 605         for (j = 0; (j < nf2); j++)
 606         {
 607             csum[j] = 0.0;
 608         }
 609         for (m = 0; (m < DIM); m++)
 610         {
 611             /* Copy the vector data in a linear array */
 612             for (j = 0; (j < nf2); j++)
 613             {
 614                 ctmp[j] = c1[DIM*j+m];
 615             }
 616             low_do_four_core(nfour, nf2, ctmp, cfour, enNorm, FALSE);
 617             for (j = 0; (j < nf2); j++)
 618             {
 619                 csum[j] += cfour[j];
 620             }
 621         }
 622     }
 623     else
 624     {
 625         gmx_fatal(FARGS, "\nUnknown mode in do_autocorr (%d)", mode);
 626     }
 627
 628     sfree(cfour);
 629     for (j = 0; (j < nf2); j++)
 630     {
 631         c1[j] = csum[j]/(real)(nframes-j);
 632     }
 633 }
 634
 635 real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose,
 636              real tbeginfit, real tendfit, real dt, real c1[], real *fit)
 637 {
 638     real        fitparm[3];
 639     real        tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate, *sig;
 640     int         i, j, jmax, nf_int;
 641     gmx_bool    bPrint;
 642
 643     bPrint = bVerbose || bDebugMode();
 644
 645     if (bPrint)
 646     {
 647         printf("COR:\n");
 648     }
 649
 650     if (tendfit <= 0)
 651     {
 652         tendfit = ncorr*dt;
 653     }
 654     nf_int = min(ncorr, (int)(tendfit/dt));
 655     sum    = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
 656
 657     /* Estimate the correlation time for better fitting */
 658     ct_estimate = 0.5*c1[0];
 659     for (i = 1; (i < ncorr) && (c1[i] > 0); i++)
 660     {
 661         ct_estimate += c1[i];
 662     }
 663     ct_estimate *= dt/c1[0];
 664
 665     if (bPrint)
 666     {
 667         printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
 668                0.0, dt*nf_int, sum);
 669         printf("COR: Relaxation times are computed as fit to an exponential:\n");
 670         printf("COR:   %s\n", longs_ffn[fitfn]);
 671         printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit, min(ncorr*dt, tendfit));
 672     }
 673
 674     tStart = 0;
 675     if (bPrint)
 676     {
 677         printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
 678                "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)",
 679                (nfp_ffn[fitfn] >= 2) ? " a2 ()" : "",
 680                (nfp_ffn[fitfn] >= 3) ? " a3 (ps)" : "");
 681     }
 682
 683     snew(sig, ncorr);
 684
 685     if (tbeginfit > 0)
 686     {
 687         jmax = 3;
 688     }
 689     else
 690     {
 691         jmax = 1;
 692     }
 693     for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++)
 694     {
 695         /* Estimate the correlation time for better fitting */
 696         c_start     = -1;
 697         ct_estimate = 0;
 698         for (i = 0; (i < ncorr) && (dt*i < tStart || c1[i] > 0); i++)
 699         {
 700             if (c_start < 0)
 701             {
 702                 if (dt*i >= tStart)
 703                 {
 704                     c_start     = c1[i];
 705                     ct_estimate = 0.5*c1[i];
 706                 }
 707             }
 708             else
 709             {
 710                 ct_estimate += c1[i];
 711             }
 712         }
 713         if (c_start > 0)
 714         {
 715             ct_estimate *= dt/c_start;
 716         }
 717         else
 718         {
 719             /* The data is strange, so we need to choose somehting */
 720             ct_estimate = tendfit;
 721         }
 722         if (debug)
 723         {
 724             fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate);
 725         }
 726
 727         if (fitfn == effnEXP3)
 728         {
 729             fitparm[0] = 0.002*ncorr*dt;
 730             fitparm[1] = 0.95;
 731             fitparm[2] = 0.2*ncorr*dt;
 732         }
 733         else
 734         {
 735             /* Good initial guess, this increases the probability of convergence */
 736             fitparm[0] = ct_estimate;
 737             fitparm[1] = 1.0;
 738             fitparm[2] = 1.0;
 739         }
 740
 741         /* Generate more or less appropriate sigma's */
 742         for (i = 0; i < ncorr; i++)
 743         {
 744             sig[i] = sqrt(ct_estimate+dt*i);
 745         }
 746
 747         nf_int    = min(ncorr, (int)((tStart+1e-4)/dt));
 748         sum       = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
 749         tail_corr = do_lmfit(ncorr, c1, sig, dt, NULL, tStart, tendfit, oenv,
 750                              bDebugMode(), fitfn, fitparm, 0);
 751         sumtot = sum+tail_corr;
 752         if (fit && ((jmax == 1) || (j == 1)))
 753         {
 754             for (i = 0; (i < ncorr); i++)
 755             {
 756                 fit[i] = fit_function(fitfn, fitparm, i*dt);
 757             }
 758         }
 759         if (bPrint)
 760         {
 761             printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot);
 762             for (i = 0; (i < nfp_ffn[fitfn]); i++)
 763             {
 764                 printf(" %11.4e", fitparm[i]);
 765             }
 766             printf("\n");
 767         }
 768         tStart += tbeginfit;
 769     }
 770     sfree(sig);
 771
 772     return sumtot;
 773 }
 774
 775 void low_do_autocorr(const char *fn, const output_env_t oenv, const char *title,
 776                      int nframes, int nitem, int nout, real **c1,
 777                      real dt, unsigned long mode, int nrestart,
 778                      gmx_bool bAver, gmx_bool bNormalize,
 779                      gmx_bool bVerbose, real tbeginfit, real tendfit,
 780                      int eFitFn)
 781 {
 782     FILE       *fp, *gp = NULL;
 783     int         i, k, nfour;
 784     real       *csum;
 785     real       *ctmp, *fit;
 786     real        c0, sum, Ct2av, Ctav;
 787     gmx_bool    bFour = acf.bFour;
 788
 789     /* Check flags and parameters */
 790     nout = get_acfnout();
 791     if (nout == -1)
 792     {
 793         nout = acf.nout = (nframes+1)/2;
 794     }
 795     else if (nout > nframes)
 796     {
 797         nout = nframes;
 798     }
 799
 800     if (MODE(eacCos) && MODE(eacVector))
 801     {
 802         gmx_fatal(FARGS, "Incompatible options bCos && bVector (%s, %d)",
 803                   __FILE__, __LINE__);
 804     }
 805     if ((MODE(eacP3) || MODE(eacRcross)) && bFour)
 806     {
 807         fprintf(stderr, "Can't combine mode %lu with FFT, turning off FFT\n", mode);
 808         bFour = FALSE;
 809     }
 810     if (MODE(eacNormal) && MODE(eacVector))
 811     {
 812         gmx_fatal(FARGS, "Incompatible mode bits: normal and vector (or Legendre)");
 813     }
 814
 815     /* Print flags and parameters */
 816     if (bVerbose)
 817     {
 818         printf("Will calculate %s of %d thingies for %d frames\n",
 819                title ? title : "autocorrelation", nitem, nframes);
 820         printf("bAver = %s, bFour = %s bNormalize= %s\n",
 821                bool_names[bAver], bool_names[bFour], bool_names[bNormalize]);
 822         printf("mode = %lu, dt = %g, nrestart = %d\n", mode, dt, nrestart);
 823     }
 824     if (bFour)
 825     {
 826         c0 = log((double)nframes)/log(2.0);
 827         k  = c0;
 828         if (k < c0)
 829         {
 830             k++;
 831         }
 832         k++;
 833         nfour = 1<<k;
 834         if (debug)
 835         {
 836             fprintf(debug, "Using FFT to calculate %s, #points for FFT = %d\n",
 837                     title, nfour);
 838         }
 839
 840         /* Allocate temp arrays */
 841         snew(csum, nfour);
 842         snew(ctmp, nfour);
 843     }
 844     else
 845     {
 846         nfour = 0; /* To keep the compiler happy */
 847         snew(csum, nframes);
 848         snew(ctmp, nframes);
 849     }
 850
 851     /* Loop over items (e.g. molecules or dihedrals)
 852      * In this loop the actual correlation functions are computed, but without
 853      * normalizing them.
 854      */
 855     k = max(1, pow(10, (int)(log(nitem)/log(100))));
 856     for (i = 0; i < nitem; i++)
 857     {
 858         if (bVerbose && ((i%k == 0 || i == nitem-1)))
 859         {
 860             fprintf(stderr, "\rThingie %d", i+1);
 861         }
 862
 863         if (bFour)
 864         {
 865             do_four_core(mode, nfour, nframes, nframes, c1[i], csum, ctmp);
 866         }
 867         else
 868         {
 869             do_ac_core(nframes, nout, ctmp, c1[i], nrestart, mode);
 870         }
 871     }
 872     if (bVerbose)
 873     {
 874         fprintf(stderr, "\n");
 875     }
 876     sfree(ctmp);
 877     sfree(csum);
 878
 879     if (fn)
 880     {
 881         snew(fit, nout);
 882         fp = xvgropen(fn, title, "Time (ps)", "C(t)", oenv);
 883     }
 884     else
 885     {
 886         fit = NULL;
 887         fp  = NULL;
 888     }
 889     if (bAver)
 890     {
 891         if (nitem > 1)
 892         {
 893             average_acf(bVerbose, nframes, nitem, c1);
 894         }
 895
 896         if (bNormalize)
 897         {
 898             normalize_acf(nout, c1[0]);
 899         }
 900
 901         if (eFitFn != effnNONE)
 902         {
 903             fit_acf(nout, eFitFn, oenv, fn != NULL, tbeginfit, tendfit, dt, c1[0], fit);
 904             sum = print_and_integrate(fp, nout, dt, c1[0], fit, 1);
 905         }
 906         else
 907         {
 908             sum = print_and_integrate(fp, nout, dt, c1[0], NULL, 1);
 909             if (bVerbose)
 910             {
 911                 printf("Correlation time (integral over corrfn): %g (ps)\n", sum);
 912             }
 913         }
 914     }
 915     else
 916     {
 917         /* Not averaging. Normalize individual ACFs */
 918         Ctav = Ct2av = 0;
 919         if (debug)
 920         {
 921             gp = xvgropen("ct-distr.xvg", "Correlation times", "item", "time (ps)", oenv);
 922         }
 923         for (i = 0; i < nitem; i++)
 924         {
 925             if (bNormalize)
 926             {
 927                 normalize_acf(nout, c1[i]);
 928             }
 929             if (eFitFn != effnNONE)
 930             {
 931                 fit_acf(nout, eFitFn, oenv, fn != NULL, tbeginfit, tendfit, dt, c1[i], fit);
 932                 sum = print_and_integrate(fp, nout, dt, c1[i], fit, 1);
 933             }
 934             else
 935             {
 936                 sum = print_and_integrate(fp, nout, dt, c1[i], NULL, 1);
 937                 if (debug)
 938                 {
 939                     fprintf(debug,
 940                             "CORRelation time (integral over corrfn %d): %g (ps)\n",
 941                             i, sum);
 942                 }
 943             }
 944             Ctav  += sum;
 945             Ct2av += sum*sum;
 946             if (debug)
 947             {
 948                 fprintf(gp, "%5d  %.3f\n", i, sum);
 949             }
 950         }
 951         if (debug)
 952         {
 953             gmx_ffclose(gp);
 954         }
 955         if (nitem > 1)
 956         {
 957             Ctav  /= nitem;
 958             Ct2av /= nitem;
 959             printf("Average correlation time %.3f Std. Dev. %.3f Error %.3f (ps)\n",
 960                    Ctav, sqrt((Ct2av - sqr(Ctav))),
 961                    sqrt((Ct2av - sqr(Ctav))/(nitem-1)));
 962         }
 963     }
 964     if (fp)
 965     {
 966         gmx_ffclose(fp);
 967     }
 968     sfree(fit);
 969 }
 970
 971 static const char *Leg[]   = { NULL, "0", "1", "2", "3", NULL };
 972
 973 t_pargs *add_acf_pargs(int *npargs, t_pargs *pa)
 974 {
 975     t_pargs  acfpa[] = {
 976         { "-acflen",     FALSE, etINT,  {&acf.nout},
 977           "Length of the ACF, default is half the number of frames" },
 978         { "-normalize", FALSE, etBOOL, {&acf.bNormalize},
 979           "Normalize ACF" },
 980         { "-fftcorr",  FALSE, etBOOL, {&acf.bFour},
 981           "HIDDENUse fast fourier transform for correlation function" },
 982         { "-nrestart", FALSE, etINT,  {&acf.nrestart},
 983           "HIDDENNumber of frames between time origins for ACF when no FFT is used" },
 984         { "-P",        FALSE, etENUM, {Leg},
 985           "Order of Legendre polynomial for ACF (0 indicates none)" },
 986         { "-fitfn",    FALSE, etENUM, {s_ffn},
 987           "Fit function" },
 988         { "-beginfit", FALSE, etREAL, {&acf.tbeginfit},
 989           "Time where to begin the exponential fit of the correlation function" },
 990         { "-endfit",   FALSE, etREAL, {&acf.tendfit},
 991           "Time where to end the exponential fit of the correlation function, -1 is until the end" },
 992     };
 993 #define NPA asize(acfpa)
 994     t_pargs *ppa;
 995     int      i;
 996
 997     snew(ppa, *npargs+NPA);
 998     for (i = 0; (i < *npargs); i++)
 999     {
1000         ppa[i] = pa[i];
1001     }
1002     for (i = 0; (i < NPA); i++)
1003     {
1004         ppa[*npargs+i] = acfpa[i];
1005     }
1006     (*npargs) += NPA;
1007
1008     acf.mode       = 0;
1009     acf.nrestart   = 1;
1010     acf.nout       = -1;
1011     acf.P          = 0;
1012     acf.fitfn      = effnEXP1;
1013     acf.bFour      = TRUE;
1014     acf.bNormalize = TRUE;
1015     acf.tbeginfit  = 0.0;
1016     acf.tendfit    = -1;
1017
1018     bACFinit = TRUE;
1019
1020     return ppa;
1021 }
1022
1023 void do_autocorr(const char *fn, const output_env_t oenv, const char *title,
1024                  int nframes, int nitem, real **c1,
1025                  real dt, unsigned long mode, gmx_bool bAver)
1026 {
1027     if (!bACFinit)
1028     {
1029         printf("ACF data structures have not been initialised. Call add_acf_pargs\n");
1030     }
1031
1032     /* Handle enumerated types */
1033     sscanf(Leg[0], "%d", &acf.P);
1034     acf.fitfn = sffn2effn(s_ffn);
1035
1036     switch (acf.P)
1037     {
1038         case 1:
1039             mode = mode | eacP1;
1040             break;
1041         case 2:
1042             mode = mode | eacP2;
1043             break;
1044         case 3:
1045             mode = mode | eacP3;
1046             break;
1047         default:
1048             break;
1049     }
1050
1051     low_do_autocorr(fn, oenv, title, nframes, nitem, acf.nout, c1, dt, mode,
1052                     acf.nrestart, bAver, acf.bNormalize,
1053                     bDebugMode(), acf.tbeginfit, acf.tendfit,
1054                     acf.fitfn);
1055 }
1056
1057 int get_acfnout(void)
1058 {
1059     if (!bACFinit)
1060     {
1061         gmx_fatal(FARGS, "ACF data not initialized yet");
1062     }
1063
1064     return acf.nout;
1065 }
1066
1067 int get_acffitfn(void)
1068 {
1069     if (!bACFinit)
1070     {
1071         gmx_fatal(FARGS, "ACF data not initialized yet");
1072     }
1073
1074     return sffn2effn(s_ffn);
1075 }
1076
1077 void cross_corr(int n, real f[], real g[], real corr[])
1078 {
1079     int i, j;
1080
1081     if (acf.bFour)
1082     {
1083         correl(f-1, g-1, n, corr-1);
1084     }
1085     else
1086     {
1087         for (i = 0; (i < n); i++)
1088         {
1089             for (j = i; (j < n); j++)
1090             {
1091                 corr[j-i] += f[i]*g[j];
1092             }
1093             corr[i] /= (real)(n-i);
1094         }
1095     }
1096 }