src/tools/autocorr.c

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