src/tools/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdio.h>
  43 #include <math.h>
  44 #include "macros.h"
  45 #include "typedefs.h"
  46 #include "physics.h"
  47 #include "smalloc.h"
  48 #include "xvgr.h"
  49 #include "futil.h"
  50 #include "gstat.h"
  51 #include "names.h"
  52 #include "gmx_fatal.h"
  53 #include "vec.h"
  54 #include "string2.h"
  55 #include "correl.h"
  56
  57 #define MODE(x) ((mode & (x)) == (x))
  58
  59 typedef struct {
  60   unsigned long mode;
  61   int  nrestart,nout,P,fitfn,nskip;
  62   gmx_bool bFour,bNormalize;
  63   real tbeginfit,tendfit;
  64 } t_acf;
  65
  66 static gmx_bool  bACFinit = FALSE;
  67 static t_acf acf;
  68
  69 enum { enNorm, enCos, enSin };
  70
  71 int sffn2effn(const char **sffn)
  72 {
  73   int eFitFn,i;
  74
  75   eFitFn = 0;
  76   for(i=0; i<effnNR; i++)
  77     if (sffn[i+1] && strcmp(sffn[0],sffn[i+1])==0)
  78       eFitFn = i;
  79
  80   return eFitFn;
  81 }
  82
  83 static void low_do_four_core(int nfour,int nframes,real c1[],real cfour[],
  84                              int nCos,gmx_bool bPadding)
  85 {
  86   int  i=0;
  87   real aver,*ans;
  88
  89   aver = 0.0;
  90   switch (nCos) {
  91   case enNorm:
  92     for(i=0; (i<nframes); i++) {
  93       aver+=c1[i];
  94       cfour[i]=c1[i];
  95     }
  96     break;
  97   case enCos:
  98     for(i=0; (i<nframes); i++)
  99       cfour[i]=cos(c1[i]);
 100     break;
 101   case enSin:
 102     for(i=0; (i<nframes); i++)
 103       cfour[i]=sin(c1[i]);
 104     break;
 105   default:
 106     gmx_fatal(FARGS,"nCos = %d, %s %d",nCos,__FILE__,__LINE__);
 107   }
 108   for(   ; (i<nfour); i++)
 109     cfour[i]= 0.0;
 110
 111   if (bPadding) {
 112     aver /= nframes;
 113     /* printf("aver = %g\n",aver); */
 114     for(i=0; (i<nframes); i++)
 115       cfour[i] -= aver;
 116   }
 117
 118   snew(ans,2*nfour);
 119   correl(cfour-1,cfour-1,nfour,ans-1);
 120
 121   if (bPadding)
 122     for (i=0; (i<nfour); i++)
 123       cfour[i] = ans[i]+sqr(aver);
 124   else
 125     for (i=0; (i<nfour); i++)
 126       cfour[i] = ans[i];
 127
 128   sfree(ans);
 129 }
 130
 131 static void do_ac_core(int nframes,int nout,
 132                        real corr[],real c1[],int nrestart,
 133                        unsigned long mode)
 134 {
 135   int     j,k,j3,jk3,m,n;
 136   real    ccc,cth;
 137   rvec    xj,xk,rr;
 138
 139   if (nrestart < 1) {
 140     printf("WARNING: setting number of restarts to 1\n");
 141     nrestart = 1;
 142   }
 143   if (debug)
 144     fprintf(debug,
 145             "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n",
 146             nframes,nout,nrestart,mode);
 147
 148   for(j=0; (j<nout); j++)
 149     corr[j]=0;
 150
 151   /* Loop over starting points. */
 152   for(j=0; (j<nframes); j+=nrestart) {
 153     j3  = DIM*j;
 154
 155     /* Loop over the correlation length for this starting point */
 156     for(k=0; (k<nout) && (j+k < nframes); k++) {
 157       jk3 = DIM*(j+k);
 158
 159       /* Switch over possible ACF types.
 160        * It might be more efficient to put the loops inside the switch,
 161        * but this is more clear, and save development time!
 162        */
 163       if (MODE(eacNormal)) {
 164         corr[k] += c1[j]*c1[j+k];
 165       }
 166       else if (MODE(eacCos)) {
 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         /* Check equality (phi(t)==phi(t+dt)) */
 172         corr[k] += (c1[j]==c1[j+k])? 1 : 0;
 173       }
 174       else if (MODE(eacP1) || MODE(eacP2) || MODE(eacP3)) {
 175         for(m=0; (m<DIM); m++) {
 176           xj[m] = c1[j3+m];
 177           xk[m] = c1[jk3+m];
 178         }
 179         cth=cos_angle(xj,xk);
 180
 181         if (cth-1.0 > 1.0e-15) {
 182           printf("j: %d, k: %d, xj:(%g,%g,%g), xk:(%g,%g,%g)\n",
 183                   j,k,xj[XX],xj[YY],xj[ZZ],xk[XX],xk[YY],xk[ZZ]);
 184         }
 185
 186         corr[k] += LegendreP(cth,mode);  /* 1.5*cth*cth-0.5; */
 187       }
 188       else if (MODE(eacRcross)) {
 189         for(m=0; (m<DIM); m++) {
 190           xj[m] = c1[j3+m];
 191           xk[m] = c1[jk3+m];
 192         }
 193         cprod(xj,xk,rr);
 194
 195         corr[k] += iprod(rr,rr);
 196       }
 197       else if (MODE(eacVector)) {
 198         for(m=0; (m<DIM); m++) {
 199           xj[m] = c1[j3+m];
 200           xk[m] = c1[jk3+m];
 201         }
 202         ccc = iprod(xj,xk);
 203
 204         corr[k] += ccc;
 205       }
 206       else
 207         gmx_fatal(FARGS,"\nInvalid mode (%d) in do_ac_core",mode);
 208     }
 209   }
 210   /* Correct for the number of points and copy results to the data array */
 211   for(j=0; (j<nout); j++) {
 212     n = (nframes-j+(nrestart-1))/nrestart;
 213     c1[j] = corr[j]/n;
 214   }
 215 }
 216
 217 void normalize_acf(int nout,real corr[])
 218 {
 219   int  j;
 220   real c0;
 221
 222   if (debug) {
 223     fprintf(debug,"Before normalization\n");
 224     for(j=0; (j<nout); j++)
 225       fprintf(debug,"%5d  %10f\n",j,corr[j]);
 226   }
 227
 228   /* Normalisation makes that c[0] = 1.0 and that other points are scaled
 229    * accordingly.
 230    */
 231   if (corr[0] == 0.0)
 232     c0 = 1.0;
 233   else
 234     c0 = 1.0/corr[0];
 235   for(j=0; (j<nout); j++)
 236     corr[j] *= c0;
 237
 238   if (debug) {
 239     fprintf(debug,"After normalization\n");
 240     for(j=0; (j<nout); j++)
 241       fprintf(debug,"%5d  %10f\n",j,corr[j]);
 242   }
 243 }
 244
 245 void average_acf(gmx_bool bVerbose,int n,int nitem,real **c1)
 246 {
 247   real c0;
 248   int  i,j;
 249
 250   if (bVerbose)
 251     printf("Averaging correlation functions\n");
 252
 253   for(j=0; (j<n); j++) {
 254     c0 = 0;
 255     for(i=0; (i<nitem); i++)
 256       c0+=c1[i][j];
 257     c1[0][j] = c0/nitem;
 258   }
 259 }
 260
 261 void norm_and_scale_vectors(int nframes,real c1[],real scale)
 262 {
 263   int  j,m;
 264   real *rij;
 265
 266   for(j=0; (j<nframes); j++) {
 267     rij = &(c1[j*DIM]);
 268     unitv(rij,rij);
 269     for(m=0; (m<DIM); m++)
 270       rij[m]*=scale;
 271   }
 272 }
 273
 274 void dump_tmp(char *s,int n,real c[])
 275 {
 276   FILE *fp;
 277   int  i;
 278
 279   fp=ffopen(s,"w");
 280   for(i=0; (i<n); i++)
 281     fprintf(fp,"%10d  %10g\n",i,c[i]);
 282   ffclose(fp);
 283 }
 284
 285 real print_and_integrate(FILE *fp,int n,real dt,real c[],real *fit,int nskip)
 286 {
 287   real c0,sum;
 288   int  j;
 289
 290   /* Use trapezoidal rule for calculating integral */
 291   sum = 0.0;
 292   for(j=0; (j<n); j++) {
 293     c0 = c[j];
 294     if (fp && (nskip == 0 || j % nskip == 0))
 295       fprintf(fp,"%10.3f  %10.5f\n",j*dt,c0);
 296     if (j > 0)
 297       sum+=dt*(c0+c[j-1]);
 298   }
 299   if (fp) {
 300     fprintf(fp,"&\n");
 301     if (fit) {
 302       for(j=0; (j<n); j++)
 303         if (nskip == 0 || j % nskip == 0)
 304           fprintf(fp,"%10.3f  %10.5f\n",j*dt,fit[j]);
 305       fprintf(fp,"&\n");
 306     }
 307   }
 308   return sum*0.5;
 309 }
 310
 311 real evaluate_integral(int n,real x[],real y[],real dy[],real aver_start,
 312                        real *stddev)
 313 {
 314   double sum,sum_var,w;
 315   double sum_tail=0,sum2_tail=0;
 316   int    j,nsum_tail=0;
 317
 318   /* Use trapezoidal rule for calculating integral */
 319   if (n <= 0)
 320     gmx_fatal(FARGS,"Evaluating integral: n = %d (file %s, line %d)",
 321               n,__FILE__,__LINE__);
 322
 323   sum = 0;
 324   sum_var = 0;
 325   for(j=0; (j<n); j++) {
 326     w = 0;
 327     if (j > 0)
 328       w += 0.5*(x[j] - x[j-1]);
 329     if (j < n-1)
 330       w += 0.5*(x[j+1] - x[j]);
 331     sum += w*y[j];
 332     if (dy) {
 333       /* Assume all errors are uncorrelated */
 334       sum_var += sqr(w*dy[j]);
 335     }
 336
 337     if ((aver_start > 0) && (x[j] >= aver_start)) {
 338       sum_tail  += sum;
 339       sum2_tail += sqrt(sum_var);
 340       nsum_tail += 1;
 341     }
 342   }
 343
 344   if (nsum_tail > 0) {
 345     sum = sum_tail/nsum_tail;
 346     /* This is a worst case estimate, assuming all stddev's are correlated. */
 347     *stddev = sum2_tail/nsum_tail;
 348   }
 349   else
 350     *stddev = sqrt(sum_var);
 351
 352   return sum;
 353 }
 354
 355 void do_four_core(unsigned long mode,int nfour,int nf2,int nframes,
 356                   real c1[],real csum[],real ctmp[])
 357 {
 358   real *cfour;
 359   char    buf[32];
 360   real    fac;
 361   int     j,m,m1;
 362
 363   snew(cfour,nfour);
 364
 365   if (MODE(eacNormal)) {
 366     /********************************************
 367      *  N O R M A L
 368      ********************************************/
 369     low_do_four_core(nfour,nf2,c1,csum,enNorm,FALSE);
 370   }
 371   else if (MODE(eacCos)) {
 372     /***************************************************
 373      * C O S I N E
 374      ***************************************************/
 375     /* Copy the data to temp array. Since we need it twice
 376      * we can't overwrite original.
 377      */
 378     for(j=0; (j<nf2); j++)
 379       ctmp[j]=c1[j];
 380
 381     /* Cosine term of AC function */
 382     low_do_four_core(nfour,nf2,ctmp,cfour,enCos,FALSE);
 383     for(j=0; (j<nf2); j++)
 384       c1[j]  = cfour[j];
 385
 386     /* Sine term of AC function */
 387     low_do_four_core(nfour,nf2,ctmp,cfour,enSin,FALSE);
 388     for(j=0; (j<nf2); j++) {
 389       c1[j] += cfour[j];
 390       csum[j] = c1[j];
 391     }
 392   }
 393   else if (MODE(eacP2)) {
 394     /***************************************************
 395      * Legendre polynomials
 396      ***************************************************/
 397     /* First normalize the vectors */
 398     norm_and_scale_vectors(nframes,c1,1.0);
 399
 400     /* For P2 thingies we have to do six FFT based correls
 401      * First for XX^2, then for YY^2, then for ZZ^2
 402      * Then we have to do XY, YZ and XZ (counting these twice)
 403      * After that we sum them and normalise
 404      * P2(x) = (3 * cos^2 (x) - 1)/2
 405      * for unit vectors u and v we compute the cosine as the inner product
 406      * cos(u,v) = uX vX + uY vY + uZ vZ
 407      *
 408      *        oo
 409      *        /
 410      * C(t) = |  (3 cos^2(u(t'),u(t'+t)) - 1)/2 dt'
 411      *        /
 412      *        0
 413      *
 414      * For ACF we need:
 415      * P2(u(0),u(t)) = [3 * (uX(0) uX(t) +
 416      *                       uY(0) uY(t) +
 417      *                       uZ(0) uZ(t))^2 - 1]/2
 418      *               = [3 * ((uX(0) uX(t))^2 +
 419      *                       (uY(0) uY(t))^2 +
 420      *                       (uZ(0) uZ(t))^2 +
 421      *                 2(uX(0) uY(0) uX(t) uY(t)) +
 422      *                 2(uX(0) uZ(0) uX(t) uZ(t)) +
 423      *                 2(uY(0) uZ(0) uY(t) uZ(t))) - 1]/2
 424      *
 425      *               = [(3/2) * (<uX^2> + <uY^2> + <uZ^2> +
 426      *                         2<uXuY> + 2<uXuZ> + 2<uYuZ>) - 0.5]
 427      *
 428      */
 429
 430     /* Because of normalization the number of -0.5 to subtract
 431      * depends on the number of data points!
 432      */
 433     for(j=0; (j<nf2); j++)
 434       csum[j]  = -0.5*(nf2-j);
 435
 436     /***** DIAGONAL ELEMENTS ************/
 437     for(m=0; (m<DIM); m++) {
 438       /* Copy the vector data in a linear array */
 439       for(j=0; (j<nf2); j++)
 440         ctmp[j]  = sqr(c1[DIM*j+m]);
 441       if (debug) {
 442         sprintf(buf,"c1diag%d.xvg",m);
 443         dump_tmp(buf,nf2,ctmp);
 444       }
 445
 446       low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
 447
 448       if (debug) {
 449         sprintf(buf,"c1dfout%d.xvg",m);
 450         dump_tmp(buf,nf2,cfour);
 451       }
 452       fac = 1.5;
 453       for(j=0; (j<nf2); j++)
 454         csum[j] += fac*(cfour[j]);
 455     }
 456     /******* OFF-DIAGONAL ELEMENTS **********/
 457     for(m=0; (m<DIM); m++) {
 458       /* Copy the vector data in a linear array */
 459       m1=(m+1) % DIM;
 460       for(j=0; (j<nf2); j++)
 461         ctmp[j]=c1[DIM*j+m]*c1[DIM*j+m1];
 462
 463       if (debug) {
 464         sprintf(buf,"c1off%d.xvg",m);
 465         dump_tmp(buf,nf2,ctmp);
 466       }
 467       low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
 468       if (debug) {
 469         sprintf(buf,"c1ofout%d.xvg",m);
 470         dump_tmp(buf,nf2,cfour);
 471       }
 472       fac = 3.0;
 473       for(j=0; (j<nf2); j++) {
 474         csum[j] += fac*cfour[j];
 475       }
 476     }
 477   }
 478   else if (MODE(eacP1) || MODE(eacVector)) {
 479     /***************************************************
 480      * V E C T O R & P1
 481      ***************************************************/
 482     if (MODE(eacP1)) {
 483       /* First normalize the vectors */
 484       norm_and_scale_vectors(nframes,c1,1.0);
 485     }
 486
 487     /* For vector thingies we have to do three FFT based correls
 488      * First for XX, then for YY, then for ZZ
 489      * After that we sum them and normalise
 490      */
 491     for(j=0; (j<nf2); j++) {
 492       csum[j]=0.0;
 493     }
 494     for(m=0; (m<DIM); m++) {
 495       /* Copy the vector data in a linear array */
 496       for(j=0; (j<nf2); j++)
 497         ctmp[j]=c1[DIM*j+m];
 498       low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
 499       for(j=0; (j<nf2); j++)
 500         csum[j] += cfour[j];
 501     }
 502   }
 503   else
 504     gmx_fatal(FARGS,"\nUnknown mode in do_autocorr (%d)",mode);
 505
 506   sfree(cfour);
 507   for(j=0; (j<nf2); j++)
 508     c1[j] = csum[j]/(real)(nframes-j);
 509 }
 510
 511 real fit_acf(int ncorr,int fitfn,const output_env_t oenv,gmx_bool bVerbose,
 512              real tbeginfit,real tendfit,real dt,real c1[],real *fit)
 513 {
 514   real    fitparm[3];
 515   real    tStart,tail_corr,sum,sumtot=0,c_start,ct_estimate,*sig;
 516   int     i,j,jmax,nf_int;
 517   gmx_bool    bPrint;
 518
 519   bPrint = bVerbose || bDebugMode();
 520
 521   if (bPrint) printf("COR:\n");
 522
 523   if (tendfit <= 0)
 524     tendfit = ncorr*dt;
 525   nf_int = min(ncorr,(int)(tendfit/dt));
 526   sum    = print_and_integrate(debug,nf_int,dt,c1,NULL,1);
 527
 528   /* Estimate the correlation time for better fitting */
 529   ct_estimate = 0.5*c1[0];
 530   for(i=1; (i<ncorr) && (c1[i]>0); i++)
 531       ct_estimate += c1[i];
 532   ct_estimate *= dt/c1[0];
 533
 534   if (bPrint) {
 535     printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
 536            0.0,dt*nf_int,sum);
 537     printf("COR: Relaxation times are computed as fit to an exponential:\n");
 538     printf("COR:   %s\n",longs_ffn[fitfn]);
 539     printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n",tbeginfit,min(ncorr*dt,tendfit));
 540   }
 541
 542   tStart = 0;
 543   if (bPrint)
 544     printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
 545            "Fit from","Integral","Tail Value","Sum (ps)"," a1 (ps)",
 546            (nfp_ffn[fitfn]>=2) ? " a2 ()" : "",
 547            (nfp_ffn[fitfn]>=3) ? " a3 (ps)" : "");
 548
 549   snew(sig,ncorr);
 550
 551   if (tbeginfit > 0) {
 552     jmax = 3;
 553   } else {
 554     jmax = 1;
 555   }
 556   for(j=0; ((j<jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++) {
 557     /* Estimate the correlation time for better fitting */
 558     c_start = -1;
 559     ct_estimate = 0;
 560     for(i=0; (i<ncorr) && (dt*i < tStart || c1[i]>0); i++) {
 561       if (c_start < 0) {
 562         if (dt*i >= tStart) {
 563           c_start     = c1[i];
 564           ct_estimate = 0.5*c1[i];
 565         }
 566       } else {
 567         ct_estimate += c1[i];
 568       }
 569     }
 570     if (c_start > 0) {
 571       ct_estimate *= dt/c_start;
 572     } else {
 573       /* The data is strange, so we need to choose somehting */
 574       ct_estimate = tendfit;
 575     }
 576     if (debug) {
 577       fprintf(debug,"tStart %g ct_estimate: %g\n",tStart,ct_estimate);
 578     }
 579
 580     if (fitfn == effnEXP3) {
 581       fitparm[0] = 0.002*ncorr*dt;
 582       fitparm[1] = 0.95;
 583       fitparm[2] = 0.2*ncorr*dt;
 584     } else {
 585       /* Good initial guess, this increases the probability of convergence */
 586       fitparm[0] = ct_estimate;
 587       fitparm[1] = 1.0;
 588       fitparm[2] = 1.0;
 589     }
 590
 591     /* Generate more or less appropriate sigma's */
 592     for(i=0; i<ncorr; i++) {
 593       sig[i] = sqrt(ct_estimate+dt*i);
 594     }
 595
 596     nf_int = min(ncorr,(int)((tStart+1e-4)/dt));
 597     sum    = print_and_integrate(debug,nf_int,dt,c1,NULL,1);
 598     tail_corr = do_lmfit(ncorr,c1,sig,dt,NULL,tStart,tendfit,oenv,
 599                          bDebugMode(),fitfn,fitparm,0);
 600     sumtot = sum+tail_corr;
 601     if (fit && ((jmax == 1) || (j == 1)))
 602       for(i=0; (i<ncorr); i++)
 603         fit[i] = fit_function(fitfn,fitparm,i*dt);
 604     if (bPrint) {
 605       printf("COR:%11.4e%11.4e%11.4e%11.4e",tStart,sum,tail_corr,sumtot);
 606       for(i=0; (i<nfp_ffn[fitfn]); i++)
 607         printf(" %11.4e",fitparm[i]);
 608       printf("\n");
 609     }
 610     tStart += tbeginfit;
 611   }
 612   sfree(sig);
 613
 614   return sumtot;
 615 }
 616
 617 void low_do_autocorr(const char *fn,const output_env_t oenv,const char *title,
 618                      int nframes,int nitem,int nout,real **c1,
 619                      real dt,unsigned long mode,int nrestart,
 620                      gmx_bool bAver,gmx_bool bNormalize,
 621                      gmx_bool bVerbose,real tbeginfit,real tendfit,
 622                      int eFitFn,int nskip)
 623 {
 624   FILE    *fp,*gp=NULL;
 625   int     i,k,nfour;
 626   real    *csum;
 627   real    *ctmp,*fit;
 628   real    c0,sum,Ct2av,Ctav;
 629   gmx_bool    bFour = acf.bFour;
 630
 631   /* Check flags and parameters */
 632   nout = get_acfnout();
 633   if (nout == -1)
 634     nout = acf.nout = (nframes+1)/2;
 635   else if (nout > nframes)
 636     nout=nframes;
 637
 638   if (MODE(eacCos) && MODE(eacVector))
 639     gmx_fatal(FARGS,"Incompatible options bCos && bVector (%s, %d)",
 640                 __FILE__,__LINE__);
 641   if ((MODE(eacP3) || MODE(eacRcross)) && bFour) {
 642     fprintf(stderr,"Can't combine mode %lu with FFT, turning off FFT\n",mode);
 643     bFour = FALSE;
 644   }
 645   if (MODE(eacNormal) && MODE(eacVector))
 646     gmx_fatal(FARGS,"Incompatible mode bits: normal and vector (or Legendre)");
 647
 648   /* Print flags and parameters */
 649   if (bVerbose) {
 650     printf("Will calculate %s of %d thingies for %d frames\n",
 651            title ? title : "autocorrelation",nitem,nframes);
 652     printf("bAver = %s, bFour = %s bNormalize= %s\n",
 653            bool_names[bAver],bool_names[bFour],bool_names[bNormalize]);
 654     printf("mode = %lu, dt = %g, nrestart = %d\n",mode,dt,nrestart);
 655   }
 656   if (bFour) {
 657     c0 = log((double)nframes)/log(2.0);
 658     k  = c0;
 659     if (k < c0)
 660       k++;
 661     k++;
 662     nfour = 1<<k;
 663     if (debug)
 664       fprintf(debug,"Using FFT to calculate %s, #points for FFT = %d\n",
 665               title,nfour);
 666
 667     /* Allocate temp arrays */
 668     snew(csum,nfour);
 669     snew(ctmp,nfour);
 670   } else {
 671     nfour = 0; /* To keep the compiler happy */
 672     snew(csum,nframes);
 673     snew(ctmp,nframes);
 674   }
 675
 676   /* Loop over items (e.g. molecules or dihedrals)
 677    * In this loop the actual correlation functions are computed, but without
 678    * normalizing them.
 679    */
 680   k = max(1,pow(10,(int)(log(nitem)/log(100))));
 681   for(i=0; i<nitem; i++) {
 682     if (bVerbose && ((i%k==0 || i==nitem-1)))
 683       fprintf(stderr,"\rThingie %d",i+1);
 684
 685     if (bFour)
 686       do_four_core(mode,nfour,nframes,nframes,c1[i],csum,ctmp);
 687     else
 688       do_ac_core(nframes,nout,ctmp,c1[i],nrestart,mode);
 689   }
 690   if (bVerbose)
 691     fprintf(stderr,"\n");
 692   sfree(ctmp);
 693   sfree(csum);
 694
 695   if (fn) {
 696     snew(fit,nout);
 697     fp=xvgropen(fn,title,"Time (ps)","C(t)",oenv);
 698   } else {
 699     fit = NULL;
 700     fp  = NULL;
 701   }
 702   if (bAver) {
 703     if (nitem > 1)
 704       average_acf(bVerbose,nframes,nitem,c1);
 705
 706     if (bNormalize)
 707       normalize_acf(nout,c1[0]);
 708
 709     if (eFitFn != effnNONE) {
 710       fit_acf(nout,eFitFn,oenv,fn!=NULL,tbeginfit,tendfit,dt,c1[0],fit);
 711       sum = print_and_integrate(fp,nout,dt,c1[0],fit,1);
 712     } else {
 713       sum = print_and_integrate(fp,nout,dt,c1[0],NULL,1);
 714       if (bVerbose)
 715         printf("Correlation time (integral over corrfn): %g (ps)\n",sum);
 716     }
 717   } else {
 718     /* Not averaging. Normalize individual ACFs */
 719     Ctav = Ct2av = 0;
 720     if (debug)
 721       gp = xvgropen("ct-distr.xvg","Correlation times","item","time (ps)",oenv);
 722     for(i=0; i<nitem; i++) {
 723       if (bNormalize)
 724         normalize_acf(nout,c1[i]);
 725       if (eFitFn != effnNONE) {
 726         fit_acf(nout,eFitFn,oenv,fn!=NULL,tbeginfit,tendfit,dt,c1[i],fit);
 727         sum = print_and_integrate(fp,nout,dt,c1[i],fit,1);
 728       } else {
 729         sum = print_and_integrate(fp,nout,dt,c1[i],NULL,1);
 730         if (debug)
 731           fprintf(debug,
 732                   "CORRelation time (integral over corrfn %d): %g (ps)\n",
 733                   i,sum);
 734       }
 735       Ctav += sum;
 736       Ct2av += sum*sum;
 737       if (debug)
 738         fprintf(gp,"%5d  %.3f\n",i,sum);
 739     }
 740     if (debug)
 741       ffclose(gp);
 742     if (nitem > 1) {
 743       Ctav  /= nitem;
 744       Ct2av /= nitem;
 745       printf("Average correlation time %.3f Std. Dev. %.3f Error %.3f (ps)\n",
 746              Ctav,sqrt((Ct2av - sqr(Ctav))),
 747              sqrt((Ct2av - sqr(Ctav))/(nitem-1)));
 748     }
 749   }
 750   if (fp)
 751     ffclose(fp);
 752   sfree(fit);
 753 }
 754
 755 static const char *Leg[]   = { NULL, "0", "1", "2", "3", NULL };
 756
 757 t_pargs *add_acf_pargs(int *npargs,t_pargs *pa)
 758 {
 759   t_pargs acfpa[] = {
 760     { "-acflen",     FALSE, etINT,  {&acf.nout},
 761       "Length of the ACF, default is half the number of frames" },
 762     { "-normalize",FALSE, etBOOL, {&acf.bNormalize},
 763       "Normalize ACF" },
 764     { "-fftcorr",  FALSE, etBOOL, {&acf.bFour},
 765       "HIDDENUse fast fourier transform for correlation function" },
 766     { "-nrestart", FALSE, etINT,  {&acf.nrestart},
 767       "HIDDENNumber of frames between time origins for ACF when no FFT is used" },
 768     { "-P",        FALSE, etENUM, {Leg},
 769       "Order of Legendre polynomial for ACF (0 indicates none)" },
 770     { "-fitfn",    FALSE, etENUM, {s_ffn},
 771       "Fit function" },
 772     { "-ncskip",   FALSE, etINT,  {&acf.nskip},
 773       "Skip this many points in the output file of correlation functions" },
 774     { "-beginfit", FALSE, etREAL, {&acf.tbeginfit},
 775       "Time where to begin the exponential fit of the correlation function" },
 776     { "-endfit",   FALSE, etREAL, {&acf.tendfit},
 777       "Time where to end the exponential fit of the correlation function, -1 is until the end" },
 778    };
 779 #define NPA asize(acfpa)
 780   t_pargs *ppa;
 781   int i;
 782
 783   snew(ppa,*npargs+NPA);
 784   for(i=0; (i<*npargs); i++)
 785     ppa[i] = pa[i];
 786   for(i=0; (i<NPA); i++)
 787     ppa[*npargs+i] = acfpa[i];
 788   (*npargs) += NPA;
 789
 790   acf.mode       = 0;
 791   acf.nrestart   = 1;
 792   acf.nout       = -1;
 793   acf.P          = 0;
 794   acf.fitfn      = effnEXP1;
 795   acf.bFour      = TRUE;
 796   acf.bNormalize = TRUE;
 797   acf.tbeginfit  = 0.0;
 798   acf.tendfit    = -1;
 799
 800   bACFinit = TRUE;
 801
 802   return ppa;
 803 }
 804
 805 void do_autocorr(const char *fn,const output_env_t oenv,const char *title,
 806                  int nframes,int nitem,real **c1,
 807                  real dt,unsigned long mode,gmx_bool bAver)
 808 {
 809   if (!bACFinit) {
 810     printf("ACF data structures have not been initialised. Call add_acf_pargs\n");
 811   }
 812
 813   /* Handle enumerated types */
 814   sscanf(Leg[0],"%d",&acf.P);
 815   acf.fitfn = sffn2effn(s_ffn);
 816
 817   switch (acf.P) {
 818   case 1:
 819     mode = mode | eacP1;
 820     break;
 821   case 2:
 822     mode = mode | eacP2;
 823     break;
 824   case 3:
 825     mode = mode | eacP3;
 826     break;
 827   default:
 828     break;
 829   }
 830
 831   low_do_autocorr(fn,oenv,title,nframes,nitem,acf.nout,c1,dt,mode,
 832                   acf.nrestart,bAver,acf.bNormalize,
 833                   bDebugMode(),acf.tbeginfit,acf.tendfit,
 834                   acf.fitfn,acf.nskip);
 835 }
 836
 837 int get_acfnout(void)
 838 {
 839   if (!bACFinit)
 840     gmx_fatal(FARGS,"ACF data not initialized yet");
 841
 842   return acf.nout;
 843 }
 844
 845 int get_acffitfn(void)
 846 {
 847   if (!bACFinit)
 848     gmx_fatal(FARGS,"ACF data not initialized yet");
 849
 850   return sffn2effn(s_ffn);
 851 }
 852
 853 void cross_corr(int n,real f[],real g[],real corr[])
 854 {
 855   int i,j;
 856
 857   if (acf.bFour)
 858     correl(f-1,g-1,n,corr-1);
 859   else {
 860     for(i=0; (i<n); i++) {
 861       for(j=i; (j<n); j++)
 862         corr[j-i] += f[i]*g[j];
 863       corr[i] /= (real)(n-i);
 864     }
 865   }
 866 }