3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
49 #include "gmx_fatal.h"
54 #define MODE(x) ((mode & (x)) == (x))
58 int nrestart,nout,P,fitfn,nskip;
59 bool bFour,bNormalize;
60 real tbeginfit,tendfit;
63 static bool bACFinit = FALSE;
66 enum { enNorm, enCos, enSin };
68 int sffn2effn(const char **sffn)
73 for(i=0; i<effnNR; i++)
74 if (sffn[i+1] && strcmp(sffn[0],sffn[i+1])==0)
80 static void low_do_four_core(int nfour,int nframes,real c1[],real cfour[],
81 int nCos,bool bPadding)
89 for(i=0; (i<nframes); i++) {
95 for(i=0; (i<nframes); i++)
99 for(i=0; (i<nframes); i++)
103 gmx_fatal(FARGS,"nCos = %d, %s %d",nCos,__FILE__,__LINE__);
105 for( ; (i<nfour); i++)
110 /* printf("aver = %g\n",aver); */
111 for(i=0; (i<nframes); i++)
116 correl(cfour-1,cfour-1,nfour,ans-1);
119 for (i=0; (i<nfour); i++)
120 cfour[i] = ans[i]+sqr(aver);
122 for (i=0; (i<nfour); i++)
128 static void do_ac_core(int nframes,int nout,
129 real corr[],real c1[],int nrestart,
137 printf("WARNING: setting number of restarts to 1\n");
142 "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n",
143 nframes,nout,nrestart,mode);
145 for(j=0; (j<nout); j++)
148 /* Loop over starting points. */
149 for(j=0; (j<nframes); j+=nrestart) {
152 /* Loop over the correlation length for this starting point */
153 for(k=0; (k<nout) && (j+k < nframes); k++) {
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!
160 if (MODE(eacNormal)) {
161 corr[k] += c1[j]*c1[j+k];
163 else if (MODE(eacCos)) {
164 /* Compute the cos (phi(t)-phi(t+dt)) */
165 corr[k] += cos(c1[j]-c1[j+k]);
167 else if (MODE(eacIden)) {
168 /* Check equality (phi(t)==phi(t+dt)) */
169 corr[k] += (c1[j]==c1[j+k])? 1 : 0;
171 else if (MODE(eacP1) || MODE(eacP2) || MODE(eacP3)) {
172 for(m=0; (m<DIM); m++) {
176 cth=cos_angle(xj,xk);
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]);
183 corr[k] += LegendreP(cth,mode); /* 1.5*cth*cth-0.5; */
185 else if (MODE(eacRcross)) {
186 for(m=0; (m<DIM); m++) {
192 corr[k] += iprod(rr,rr);
194 else if (MODE(eacVector)) {
195 for(m=0; (m<DIM); m++) {
204 gmx_fatal(FARGS,"\nInvalid mode (%d) in do_ac_core",mode);
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;
214 void normalize_acf(int nout,real corr[])
220 fprintf(debug,"Before normalization\n");
221 for(j=0; (j<nout); j++)
222 fprintf(debug,"%5d %10f\n",j,corr[j]);
225 /* Normalisation makes that c[0] = 1.0 and that other points are scaled
232 for(j=0; (j<nout); j++)
236 fprintf(debug,"After normalization\n");
237 for(j=0; (j<nout); j++)
238 fprintf(debug,"%5d %10f\n",j,corr[j]);
242 void average_acf(bool bVerbose,int n,int nitem,real **c1)
248 printf("Averaging correlation functions\n");
250 for(j=0; (j<n); j++) {
252 for(i=0; (i<nitem); i++)
258 void norm_and_scale_vectors(int nframes,real c1[],real scale)
263 for(j=0; (j<nframes); j++) {
266 for(m=0; (m<DIM); m++)
271 void dump_tmp(char *s,int n,real c[])
278 fprintf(fp,"%10d %10g\n",i,c[i]);
282 real print_and_integrate(FILE *fp,int n,real dt,real c[],real *fit,int nskip)
287 /* Use trapezoidal rule for calculating integral */
289 for(j=0; (j<n); j++) {
291 if (fp && (nskip == 0 || j % nskip == 0))
292 fprintf(fp,"%10.3f %10.5f\n",j*dt,c0);
300 if (nskip == 0 || j % nskip == 0)
301 fprintf(fp,"%10.3f %10.5f\n",j*dt,fit[j]);
308 real evaluate_integral(int n,real x[],real y[],real dy[],real aver_start,
311 double sum,sum_var,w;
312 double sum_tail=0,sum2_tail=0;
315 /* Use trapezoidal rule for calculating integral */
317 gmx_fatal(FARGS,"Evaluating integral: n = %d (file %s, line %d)",
318 n,__FILE__,__LINE__);
322 for(j=0; (j<n); j++) {
325 w += 0.5*(x[j] - x[j-1]);
327 w += 0.5*(x[j+1] - x[j]);
330 /* Assume all errors are uncorrelated */
331 sum_var += sqr(w*dy[j]);
334 if ((aver_start > 0) && (x[j] >= aver_start)) {
336 sum2_tail += sqrt(sum_var);
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;
347 *stddev = sqrt(sum_var);
352 void do_four_core(unsigned long mode,int nfour,int nf2,int nframes,
353 real c1[],real csum[],real ctmp[])
362 if (MODE(eacNormal)) {
363 /********************************************
365 ********************************************/
366 low_do_four_core(nfour,nf2,c1,csum,enNorm,FALSE);
368 else if (MODE(eacCos)) {
369 /***************************************************
371 ***************************************************/
372 /* Copy the data to temp array. Since we need it twice
373 * we can't overwrite original.
375 for(j=0; (j<nf2); j++)
378 /* Cosine term of AC function */
379 low_do_four_core(nfour,nf2,ctmp,cfour,enCos,FALSE);
380 for(j=0; (j<nf2); j++)
383 /* Sine term of AC function */
384 low_do_four_core(nfour,nf2,ctmp,cfour,enSin,FALSE);
385 for(j=0; (j<nf2); j++) {
390 else if (MODE(eacP2)) {
391 /***************************************************
392 * Legendre polynomials
393 ***************************************************/
394 /* First normalize the vectors */
395 norm_and_scale_vectors(nframes,c1,1.0);
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
407 * C(t) = | (3 cos^2(u(t'),u(t'+t)) - 1)/2 dt'
412 * P2(u(0),u(t)) = [3 * (uX(0) uX(t) +
414 * uZ(0) uZ(t))^2 - 1]/2
415 * = [3 * ((uX(0) uX(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
422 * = [(3/2) * (<uX^2> + <uY^2> + <uZ^2> +
423 * 2<uXuY> + 2<uXuZ> + 2<uYuZ>) - 0.5]
427 /* Because of normalization the number of -0.5 to subtract
428 * depends on the number of data points!
430 for(j=0; (j<nf2); j++)
431 csum[j] = -0.5*(nf2-j);
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]);
439 sprintf(buf,"c1diag%d.xvg",m);
440 dump_tmp(buf,nf2,ctmp);
443 low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
446 sprintf(buf,"c1dfout%d.xvg",m);
447 dump_tmp(buf,nf2,cfour);
450 for(j=0; (j<nf2); j++)
451 csum[j] += fac*(cfour[j]);
453 /******* OFF-DIAGONAL ELEMENTS **********/
454 for(m=0; (m<DIM); m++) {
455 /* Copy the vector data in a linear array */
457 for(j=0; (j<nf2); j++)
458 ctmp[j]=c1[DIM*j+m]*c1[DIM*j+m1];
461 sprintf(buf,"c1off%d.xvg",m);
462 dump_tmp(buf,nf2,ctmp);
464 low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
466 sprintf(buf,"c1ofout%d.xvg",m);
467 dump_tmp(buf,nf2,cfour);
470 for(j=0; (j<nf2); j++) {
471 csum[j] += fac*cfour[j];
475 else if (MODE(eacP1) || MODE(eacVector)) {
476 /***************************************************
478 ***************************************************/
480 /* First normalize the vectors */
481 norm_and_scale_vectors(nframes,c1,1.0);
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
488 for(j=0; (j<nf2); j++) {
491 for(m=0; (m<DIM); m++) {
492 /* Copy the vector data in a linear array */
493 for(j=0; (j<nf2); j++)
495 low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
496 for(j=0; (j<nf2); j++)
501 gmx_fatal(FARGS,"\nUnknown mode in do_autocorr (%d)",mode);
504 for(j=0; (j<nf2); j++)
505 c1[j] = csum[j]/(real)(nframes-j);
508 real fit_acf(int ncorr,int fitfn,const output_env_t oenv,bool bVerbose,
509 real tbeginfit,real tendfit,real dt,real c1[],real *fit)
512 real tStart,tail_corr,sum,sumtot=0,ct_estimate,*sig;
516 bPrint = bVerbose || bDebugMode();
518 if (bPrint) printf("COR:\n");
522 nf_int = min(ncorr,(int)(tendfit/dt));
523 sum = print_and_integrate(debug,nf_int,dt,c1,NULL,1);
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];
532 printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
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));
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)" : "");
549 if (fitfn == effnEXP3) {
550 fitparm[0] = 0.002*ncorr*dt;
552 fitparm[2] = 0.2*ncorr*dt;
554 /* Good initial guess, this increases the probability of convergence */
555 fitparm[0] = ct_estimate;
560 /* Generate more or less appropriate sigma's */
562 for(i=0; i<ncorr; i++)
563 sig[i] = sqrt(ct_estimate+dt*i);
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);
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]);
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 bool bAver,bool bNormalize,
592 bool bVerbose,real tbeginfit,real tendfit,
593 int eFitFn,int nskip)
599 real c0,sum,Ct2av,Ctav;
600 bool bFour = acf.bFour;
602 /* Check flags and parameters */
603 nout = get_acfnout();
605 nout = acf.nout = (nframes+1)/2;
606 else if (nout > nframes)
609 if (MODE(eacCos) && MODE(eacVector))
610 gmx_fatal(FARGS,"Incompatible options bCos && bVector (%s, %d)",
612 if ((MODE(eacP3) || MODE(eacRcross)) && bFour) {
613 fprintf(stderr,"Can't combine mode %lu with FFT, turning off FFT\n",mode);
616 if (MODE(eacNormal) && MODE(eacVector))
617 gmx_fatal(FARGS,"Incompatible mode bits: normal and vector (or Legendre)");
619 /* Print flags and parameters */
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);
628 c0 = log((double)nframes)/log(2.0);
635 fprintf(debug,"Using FFT to calculate %s, #points for FFT = %d\n",
638 /* Allocate temp arrays */
642 nfour = 0; /* To keep the compiler happy */
647 /* Loop over items (e.g. molecules or dihedrals)
648 * In this loop the actual correlation functions are computed, but without
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);
657 do_four_core(mode,nfour,nframes,nframes,c1[i],csum,ctmp);
659 do_ac_core(nframes,nout,ctmp,c1[i],nrestart,mode);
662 fprintf(stderr,"\n");
668 fp=xvgropen(fn,title,"Time (ps)","C(t)",oenv);
675 average_acf(bVerbose,nframes,nitem,c1);
678 normalize_acf(nout,c1[0]);
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);
684 sum = print_and_integrate(fp,nout,dt,c1[0],NULL,1);
686 printf("Correlation time (integral over corrfn): %g (ps)\n",sum);
689 /* Not averaging. Normalize individual ACFs */
692 gp = xvgropen("ct-distr.xvg","Correlation times","item","time (ps)",oenv);
693 for(i=0; i<nitem; i++) {
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);
700 sum = print_and_integrate(fp,nout,dt,c1[i],NULL,1);
703 "CORRelation time (integral over corrfn %d): %g (ps)\n",
709 fprintf(gp,"%5d %.3f\n",i,sum);
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)));
726 static const char *Leg[] = { NULL, "0", "1", "2", "3", NULL };
728 t_pargs *add_acf_pargs(int *npargs,t_pargs *pa)
731 { "-acflen", FALSE, etINT, {&acf.nout},
732 "Length of the ACF, default is half the number of frames" },
733 { "-normalize",FALSE, etBOOL, {&acf.bNormalize},
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},
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" },
750 #define NPA asize(acfpa)
754 snew(ppa,*npargs+NPA);
755 for(i=0; (i<*npargs); i++)
757 for(i=0; (i<NPA); i++)
758 ppa[*npargs+i] = acfpa[i];
765 acf.fitfn = effnEXP1;
767 acf.bNormalize = TRUE;
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,bool bAver)
781 printf("ACF data structures have not been initialised. Call add_acf_pargs\n");
784 /* Handle enumerated types */
785 sscanf(Leg[0],"%d",&acf.P);
786 acf.fitfn = sffn2effn(s_ffn);
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);
808 int get_acfnout(void)
811 gmx_fatal(FARGS,"ACF data not initialized yet");
816 int get_acffitfn(void)
819 gmx_fatal(FARGS,"ACF data not initialized yet");
821 return sffn2effn(s_ffn);
824 void cross_corr(int n,real f[],real g[],real corr[])
829 correl(f-1,g-1,n,corr-1);
831 for(i=0; (i<n); i++) {
833 corr[j-i] += f[i]*g[j];
834 corr[i] /= (real)(n-i);