1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
46 #include "gmx_fatal.h"
54 #include "gmx_statistics.h"
58 /* must correspond to char *avbar_opt[] declared in main() */
59 enum { avbarSEL
, avbarNONE
, avbarSTDDEV
, avbarERROR
, avbar90
, avbarNR
};
61 static void power_fit(int n
,int nset
,real
**val
,real
*t
)
63 real
*x
,*y
,quality
,a
,b
,r
;
74 fprintf(stdout
,"First time is not larger than 0, using index number as time for power fit\n");
79 for(s
=0; s
<nset
; s
++) {
81 for(i
=0; i
<n
&& val
[s
][i
]>=0; i
++)
82 y
[i
] = log(val
[s
][i
]);
84 fprintf(stdout
,"Will power fit up to point %d, since it is not larger than 0\n",i
);
85 lsq_y_ax_b(i
,x
,y
,&a
,&b
,&r
,&quality
);
86 fprintf(stdout
,"Power fit set %3d: error %.3f a %g b %g\n",
87 s
+1,quality
,a
,exp(b
));
94 static real
cosine_content(int nhp
,int n
,real
*y
)
95 /* Assumes n equidistant points */
97 double fac
,cosyint
,yyint
;
103 fac
= M_PI
*nhp
/(n
-1);
108 cosyint
+= cos(fac
*i
)*y
[i
];
112 return 2*cosyint
*cosyint
/(n
*yyint
);
115 static void plot_coscont(const char *ccfile
,int n
,int nset
,real
**val
,
116 const output_env_t oenv
)
122 fp
= xvgropen(ccfile
,"Cosine content","set / half periods","cosine content",
125 for(s
=0; s
<nset
; s
++) {
126 cc
= cosine_content(s
+1,n
,val
[s
]);
127 fprintf(fp
," %d %g\n",s
+1,cc
);
128 fprintf(stdout
,"Cosine content of set %d with %.1f periods: %g\n",
131 fprintf(stdout
,"\n");
136 static void regression_analysis(int n
,bool bXYdy
,real
*x
,real
**val
)
138 real S
,chi2
,a
,b
,da
,db
,r
=0;
140 printf("Fitting data to a function f(x) = ax + b\n");
141 printf("Minimizing residual chi2 = Sum_i w_i [f(x_i) - y_i]2\n");
142 printf("Error estimates will be given if w_i (sigma) values are given\n");
143 printf("(use option -xydy).\n\n");
145 lsq_y_ax_b_error(n
,x
,val
[0],val
[1],&a
,&b
,&da
,&db
,&r
,&S
);
147 lsq_y_ax_b(n
,x
,val
[0],&a
,&b
,&r
,&S
);
149 printf("Chi2 = %g\n",chi2
);
150 printf("S (Sqrt(Chi2/(n-2)) = %g\n",S
);
151 printf("Correlation coefficient = %.1f%%\n",100*r
);
154 printf("a = %g +/- %g\n",a
,da
);
155 printf("b = %g +/- %g\n",b
,db
);
158 printf("a = %g\n",a
);
159 printf("b = %g\n",b
);
163 void histogram(const char *distfile
,real binwidth
,int n
, int nset
, real
**val
,
164 const output_env_t oenv
)
170 #if (defined SIZEOF_LONG_LONG_INT) && (SIZEOF_LONG_LONG_INT >= 8)
171 long long int *histo
;
178 for(s
=0; s
<nset
; s
++)
182 else if (val
[s
][i
] > max
)
185 min
= binwidth
*floor(min
/binwidth
);
186 max
= binwidth
*ceil(max
/binwidth
);
191 nbin
= (int)((max
- min
)/binwidth
+ 0.5) + 1;
192 fprintf(stderr
,"Making distributions with %d bins\n",nbin
);
194 fp
= xvgropen(distfile
,"Distribution","","",oenv
);
195 for(s
=0; s
<nset
; s
++) {
196 for(i
=0; i
<nbin
; i
++)
199 histo
[(int)((val
[s
][i
] - min
)/binwidth
+ 0.5)]++;
200 for(i
=0; i
<nbin
; i
++)
201 fprintf(fp
," %g %g\n",min
+i
*binwidth
,(double)histo
[i
]/(n
*binwidth
));
208 static int real_comp(const void *a
,const void *b
)
210 real dif
= *(real
*)a
- *(real
*)b
;
220 static void average(const char *avfile
,int avbar_opt
,
221 int n
, int nset
,real
**val
,real
*t
)
228 fp
= ffopen(avfile
,"w");
229 if ((avbar_opt
== avbarERROR
) && (nset
== 1))
230 avbar_opt
= avbarNONE
;
231 if (avbar_opt
!= avbarNONE
) {
232 if (avbar_opt
== avbar90
) {
234 fprintf(fp
,"@TYPE xydydy\n");
235 edge
= (int)(nset
*0.05+0.5);
236 fprintf(stdout
,"Errorbars: discarding %d points on both sides: %d%%"
237 " interval\n",edge
,(int)(100*(nset
-2*edge
)/nset
+0.5));
239 fprintf(fp
,"@TYPE xydy\n");
244 for(s
=0; s
<nset
; s
++)
247 fprintf(fp
," %g %g",t
[i
],av
);
249 if (avbar_opt
!= avbarNONE
) {
250 if (avbar_opt
== avbar90
) {
251 for(s
=0; s
<nset
; s
++)
253 qsort(tmp
,nset
,sizeof(tmp
[0]),real_comp
);
254 fprintf(fp
," %g %g",tmp
[nset
-1-edge
]-av
,av
-tmp
[edge
]);
256 for(s
=0; s
<nset
; s
++)
257 var
+= sqr(val
[s
][i
]-av
);
258 if (avbar_opt
== avbarSTDDEV
)
259 err
= sqrt(var
/nset
);
261 err
= sqrt(var
/(nset
*(nset
-1)));
262 fprintf(fp
," %g",err
);
269 if (avbar_opt
== avbar90
)
273 static real
anal_ee_inf(real
*parm
,real T
)
275 return sqrt(parm
[1]*2*parm
[0]/T
+parm
[3]*2*parm
[2]/T
);
278 static real
anal_ee(real
*parm
,real T
,real t
)
283 e1
= exp(-t
/parm
[0]);
287 e2
= exp(-t
/parm
[2]);
291 return sqrt(parm
[1]*2*parm
[0]/T
*((e1
- 1)*parm
[0]/t
+ 1) +
292 parm
[3]*2*parm
[2]/T
*((e2
- 1)*parm
[2]/t
+ 1));
295 static void estimate_error(const char *eefile
,int nb_min
,int resol
,int n
,
296 int nset
, double *av
,double *sig
,real
**val
,real dt
,
297 bool bFitAc
,bool bSingleExpFit
,bool bAllowNegLTCorr
,
298 const output_env_t oenv
)
301 int bs
,prev_bs
,nbs
,nb
;
306 real
*tbs
,*ybs
,rtmp
,dens
,*fitsig
,twooe
,tau1_est
,tau_sig
;
312 fprintf(stdout
,"The number of points is smaller than 4, can not make an error estimate\n");
317 fp
= xvgropen(eefile
,"Error estimates",
318 "Block size (time)","Error estimate", oenv
);
319 if (output_env_get_print_xvgr_codes(oenv
))
322 "@ subtitle \"using block averaging, total time %g (%d points)\"\n",
326 xvgr_legend(fp
,2*nset
,leg
,oenv
);
329 spacing
= pow(2,1.0/resol
);
333 for(s
=0; s
<nset
; s
++)
350 blav
+= val
[s
][bs
*i
+j
];
352 var
+= sqr(av
[s
] - blav
/bs
);
361 ybs
[nbs
] = var
/(nb
*(nb
-1.0))*(n
*dt
)/(sig
[s
]*sig
[s
]);
378 for(i
=0; i
<nbs
/2; i
++)
381 tbs
[i
] = tbs
[nbs
-1-i
];
384 ybs
[i
] = ybs
[nbs
-1-i
];
387 /* The initial slope of the normalized ybs^2 is 1.
388 * For a single exponential autocorrelation: ybs(tau1) = 2/e tau1
389 * From this we take our initial guess for tau1.
397 } while (i
< nbs
- 1 &&
398 (ybs
[i
] > ybs
[i
+1] || ybs
[i
] > twooe
*tau1_est
));
402 fprintf(stdout
,"Data set %d has strange time correlations:\n"
403 "the std. error using single points is larger than that of blocks of 2 points\n"
404 "The error estimate might be inaccurate, check the fit\n",
406 /* Use the total time as tau for the fitting weights */
407 tau_sig
= (n
- 1)*dt
;
416 fprintf(debug
,"set %d tau1 estimate %f\n",s
+1,tau1_est
);
419 /* Generate more or less appropriate sigma's,
420 * also taking the density of points into account.
426 dens
= tbs
[1]/tbs
[0] - 1;
430 dens
= tbs
[nbs
-1]/tbs
[nbs
-2] - 1;
434 dens
= 0.5*(tbs
[i
+1]/tbs
[i
-1] - 1);
436 fitsig
[i
] = sqrt((tau_sig
+ tbs
[i
])/dens
);
441 fitparm
[0] = tau1_est
;
443 /* We set the initial guess for tau2
444 * to halfway between tau1_est and the total time (on log scale).
446 fitparm
[2] = sqrt(tau1_est
*(n
-1)*dt
);
447 do_lmfit(nbs
,ybs
,fitsig
,0,tbs
,0,dt
*n
,oenv
,
448 bDebugMode(),effnERREST
,fitparm
,0);
449 fitparm
[3] = 1-fitparm
[1];
451 if (bSingleExpFit
|| fitparm
[0]<0 || fitparm
[2]<0 || fitparm
[1]<0
452 || (fitparm
[1]>1 && !bAllowNegLTCorr
) || fitparm
[2]>(n
-1)*dt
)
456 if (fitparm
[2] > (n
-1)*dt
)
459 "Warning: tau2 is longer than the length of the data (%g)\n"
460 " the statistics might be bad\n",
465 fprintf(stdout
,"a fitted parameter is negative\n");
467 fprintf(stdout
,"invalid fit: e.e. %g a %g tau1 %g tau2 %g\n",
468 sig
[s
]*anal_ee_inf(fitparm
,n
*dt
),
469 fitparm
[1],fitparm
[0],fitparm
[2]);
470 /* Do a fit with tau2 fixed at the total time.
471 * One could also choose any other large value for tau2.
473 fitparm
[0] = tau1_est
;
475 fitparm
[2] = (n
-1)*dt
;
476 fprintf(stderr
,"Will fix tau2 at the total time: %g\n",fitparm
[2]);
477 do_lmfit(nbs
,ybs
,fitsig
,0,tbs
,0,dt
*n
,oenv
,bDebugMode(),
478 effnERREST
,fitparm
,4);
479 fitparm
[3] = 1-fitparm
[1];
481 if (bSingleExpFit
|| fitparm
[0]<0 || fitparm
[1]<0
482 || (fitparm
[1]>1 && !bAllowNegLTCorr
))
484 if (!bSingleExpFit
) {
485 fprintf(stdout
,"a fitted parameter is negative\n");
486 fprintf(stdout
,"invalid fit: e.e. %g a %g tau1 %g tau2 %g\n",
487 sig
[s
]*anal_ee_inf(fitparm
,n
*dt
),
488 fitparm
[1],fitparm
[0],fitparm
[2]);
490 /* Do a single exponential fit */
491 fprintf(stderr
,"Will use a single exponential fit for set %d\n",s
+1);
492 fitparm
[0] = tau1_est
;
495 do_lmfit(nbs
,ybs
,fitsig
,0,tbs
,0,dt
*n
,oenv
,bDebugMode(),
496 effnERREST
,fitparm
,6);
497 fitparm
[3] = 1-fitparm
[1];
500 ee
= sig
[s
]*anal_ee_inf(fitparm
,n
*dt
);
505 fprintf(stdout
,"Set %3d: err.est. %g a %g tau1 %g tau2 %g\n",
507 fprintf(fp
,"@ legend string %d \"av %f\"\n",2*s
,av
[s
]);
508 fprintf(fp
,"@ legend string %d \"ee %6g\"\n",
509 2*s
+1,sig
[s
]*anal_ee_inf(fitparm
,n
*dt
));
512 fprintf(fp
,"%g %g %g\n",tbs
[i
],sig
[s
]*sqrt(ybs
[i
]/(n
*dt
)),
513 sig
[s
]*sqrt(fit_function(effnERREST
,fitparm
,tbs
[i
])/(n
*dt
)));
519 real
*ac
,acint
,ac_fit
[4];
523 ac
[i
] = val
[s
][i
] - av
[s
];
529 low_do_autocorr(NULL
,oenv
,NULL
,n
,1,-1,&ac
,
530 dt
,eacNormal
,1,FALSE
,TRUE
,
531 FALSE
,0,0,effnNONE
,0);
535 /* Integrate ACF only up to fitlen/2 to avoid integrating noise */
537 for(i
=1; i
<=fitlen
/2; i
++)
543 /* Generate more or less appropriate sigma's */
544 for(i
=0; i
<=fitlen
; i
++)
546 fitsig
[i
] = sqrt(acint
+ dt
*i
);
549 ac_fit
[0] = 0.5*acint
;
551 ac_fit
[2] = 10*acint
;
552 do_lmfit(n
/nb_min
,ac
,fitsig
,dt
,0,0,fitlen
*dt
,oenv
,
553 bDebugMode(),effnEXP3
,ac_fit
,0);
554 ac_fit
[3] = 1 - ac_fit
[1];
556 fprintf(stdout
,"Set %3d: ac erest %g a %g tau1 %g tau2 %g\n",
557 s
+1,sig
[s
]*anal_ee_inf(ac_fit
,n
*dt
),
558 ac_fit
[1],ac_fit
[0],ac_fit
[2]);
563 fprintf(fp
,"%g %g\n",tbs
[i
],
564 sig
[s
]*sqrt(fit_function(effnERREST
,ac_fit
,tbs
[i
]))/(n
*dt
));
580 static void luzar_correl(int nn
,real
*time
,int nset
,real
**val
,real temp
,
581 bool bError
,real fit_start
,real smooth_tail_start
,
582 const output_env_t oenv
)
584 const real tol
= 1e-8;
589 please_cite(stdout
,"Spoel2006b");
591 /* Compute negative derivative k(t) = -dc(t)/dt */
594 compute_derivative(nn
,time
,val
[0],kt
);
595 for(j
=0; (j
<nn
); j
++)
599 for(j
=0; (j
<nn
); j
++)
600 d2
+= sqr(kt
[j
] - val
[3][j
]);
601 fprintf(debug
,"RMS difference in derivatives is %g\n",sqrt(d2
/nn
));
603 analyse_corr(nn
,time
,val
[0],val
[2],kt
,NULL
,NULL
,NULL
,fit_start
,
604 temp
,smooth_tail_start
,oenv
);
607 else if (nset
== 6) {
608 analyse_corr(nn
,time
,val
[0],val
[2],val
[4],
609 val
[1],val
[3],val
[5],fit_start
,temp
,smooth_tail_start
,oenv
);
612 printf("Inconsistent input. I need c(t) sigma_c(t) n(t) sigma_n(t) K(t) sigma_K(t)\n");
613 printf("Not doing anything. Sorry.\n");
617 static void filter(real flen
,int n
,int nset
,real
**val
,real dt
,
618 const output_env_t oenv
)
621 double *filt
,sum
,vf
,fluc
,fluctot
;
623 f
= (int)(flen
/(2*dt
));
627 for(i
=1; i
<=f
; i
++) {
628 filt
[i
] = cos(M_PI
*dt
*i
/flen
);
633 fprintf(stdout
,"Will calculate the fluctuation over %d points\n",n
-2*f
);
634 fprintf(stdout
," using a filter of length %g of %d points\n",flen
,2*f
+1);
636 for(s
=0; s
<nset
; s
++) {
638 for(i
=f
; i
<n
-f
; i
++) {
639 vf
= filt
[0]*val
[s
][i
];
641 vf
+= filt
[j
]*(val
[s
][i
-f
]+val
[s
][i
+f
]);
642 fluc
+= sqr(val
[s
][i
] - vf
);
646 fprintf(stdout
,"Set %3d filtered fluctuation: %12.6e\n",s
+1,sqrt(fluc
));
648 fprintf(stdout
,"Overall filtered fluctuation: %12.6e\n",sqrt(fluctot
/nset
));
649 fprintf(stdout
,"\n");
654 static void do_fit(FILE *out
,int n
,bool bYdy
,int ny
,real
*x0
,real
**val
,
655 int npargs
,t_pargs
*ppa
,const output_env_t oenv
)
657 real
*c1
=NULL
,*sig
=NULL
,*fitparm
;
658 real dt
=0,tendfit
,tbeginfit
;
661 efitfn
= get_acffitfn();
662 nparm
= nfp_ffn
[efitfn
];
663 fprintf(out
,"Will fit to the following function:\n");
664 fprintf(out
,"%s\n",longs_ffn
[efitfn
]);
669 fprintf(out
,"Using two columns as y and sigma values\n");
673 if (opt2parg_bSet("-beginfit",npargs
,ppa
)) {
674 tbeginfit
= opt2parg_real("-beginfit",npargs
,ppa
);
676 tbeginfit
= x0
? x0
[0] : 0;
678 if (opt2parg_bSet("-endfit",npargs
,ppa
)) {
679 tendfit
= opt2parg_real("-endfit",npargs
,ppa
);
681 tendfit
= x0
? x0
[ny
-1] : (ny
-1)*dt
;
695 fitparm
[1] = 0.5*c1
[0];
699 fitparm
[0] = fitparm
[2] = 0.5*c1
[0];
705 fitparm
[0] = fitparm
[2] = fitparm
[4] = 0.33*c1
[0];
712 fitparm
[0] = fitparm
[2] = fitparm
[4] = fitparm
[6] = 0.25*c1
[0];
720 fprintf(out
,"Warning: don't know how to initialize the parameters\n");
721 for(i
=0; (i
<nparm
); i
++)
724 fprintf(out
,"Starting parameters:\n");
725 for(i
=0; (i
<nparm
); i
++)
726 fprintf(out
,"a%-2d = %12.5e\n",i
+1,fitparm
[i
]);
727 if (do_lmfit(ny
,c1
,sig
,dt
,x0
,tbeginfit
,tendfit
,
728 oenv
,bDebugMode(),efitfn
,fitparm
,0)) {
729 for(i
=0; (i
<nparm
); i
++)
730 fprintf(out
,"a%-2d = %12.5e\n",i
+1,fitparm
[i
]);
733 fprintf(out
,"No solution was found\n");
737 int gmx_analyze(int argc
,char *argv
[])
739 static const char *desc
[] = {
740 "g_analyze reads an ascii file and analyzes data sets.",
741 "A line in the input file may start with a time",
742 "(see option [TT]-time[tt]) and any number of y values may follow.",
743 "Multiple sets can also be",
744 "read when they are seperated by & (option [TT]-n[tt]),",
745 "in this case only one y value is read from each line.",
746 "All lines starting with # and @ are skipped.",
747 "All analyses can also be done for the derivative of a set",
748 "(option [TT]-d[tt]).[PAR]",
750 "All options, except for [TT]-av[tt] and [TT]-power[tt] assume that the",
751 "points are equidistant in time.[PAR]",
753 "g_analyze always shows the average and standard deviation of each",
754 "set. For each set it also shows the relative deviation of the third",
755 "and forth cumulant from those of a Gaussian distribution with the same",
756 "standard deviation.[PAR]",
758 "Option [TT]-ac[tt] produces the autocorrelation function(s).[PAR]",
760 "Option [TT]-cc[tt] plots the resemblance of set i with a cosine of",
761 "i/2 periods. The formula is:[BR]"
762 "2 (int0-T y(t) cos(i pi t) dt)^2 / int0-T y(t) y(t) dt[BR]",
763 "This is useful for principal components obtained from covariance",
764 "analysis, since the principal components of random diffusion are",
765 "pure cosines.[PAR]",
767 "Option [TT]-msd[tt] produces the mean square displacement(s).[PAR]",
769 "Option [TT]-dist[tt] produces distribution plot(s).[PAR]",
771 "Option [TT]-av[tt] produces the average over the sets.",
772 "Error bars can be added with the option [TT]-errbar[tt].",
773 "The errorbars can represent the standard deviation, the error",
774 "(assuming the points are independent) or the interval containing",
775 "90% of the points, by discarding 5% of the points at the top and",
778 "Option [TT]-ee[tt] produces error estimates using block averaging.",
779 "A set is divided in a number of blocks and averages are calculated for",
780 "each block. The error for the total average is calculated from",
781 "the variance between averages of the m blocks B_i as follows:",
782 "error^2 = Sum (B_i - <B>)^2 / (m*(m-1)).",
783 "These errors are plotted as a function of the block size.",
784 "Also an analytical block average curve is plotted, assuming",
785 "that the autocorrelation is a sum of two exponentials.",
786 "The analytical curve for the block average is:[BR]",
787 "f(t) = sigma sqrt(2/T ( a (tau1 ((exp(-t/tau1) - 1) tau1/t + 1)) +[BR]",
788 " (1-a) (tau2 ((exp(-t/tau2) - 1) tau2/t + 1)))),[BR]"
789 "where T is the total time.",
790 "a, tau1 and tau2 are obtained by fitting f^2(t) to error^2.",
791 "When the actual block average is very close to the analytical curve,",
792 "the error is sigma*sqrt(2/T (a tau1 + (1-a) tau2)).",
793 "The complete derivation is given in",
794 "B. Hess, J. Chem. Phys. 116:209-217, 2002.[PAR]",
796 "Option [TT]-filter[tt] prints the RMS high-frequency fluctuation",
797 "of each set and over all sets with respect to a filtered average.",
798 "The filter is proportional to cos(pi t/len) where t goes from -len/2",
799 "to len/2. len is supplied with the option [TT]-filter[tt].",
800 "This filter reduces oscillations with period len/2 and len by a factor",
801 "of 0.79 and 0.33 respectively.[PAR]",
803 "Option [TT]-g[tt] fits the data to the function given with option",
804 "[TT]-fitfn[tt].[PAR]",
806 "Option [TT]-power[tt] fits the data to b t^a, which is accomplished",
807 "by fitting to a t + b on log-log scale. All points after the first",
808 "zero or negative value are ignored.[PAR]"
810 "Option [TT]-luzar[tt] performs a Luzar & Chandler kinetics analysis",
811 "on output from [TT]g_hbond[tt]. The input file can be taken directly",
812 "from [TT]g_hbond -ac[tt], and then the same result should be produced."
814 static real tb
=-1,te
=-1,frac
=0.5,filtlen
=0,binwidth
=0.1,aver_start
=0;
815 static bool bHaveT
=TRUE
,bDer
=FALSE
,bSubAv
=TRUE
,bAverCorr
=FALSE
,bXYdy
=FALSE
;
816 static bool bEESEF
=FALSE
,bEENLC
=FALSE
,bEeFitAc
=FALSE
,bPower
=FALSE
;
817 static bool bIntegrate
=FALSE
,bRegression
=FALSE
,bLuzar
=FALSE
,bLuzarError
=FALSE
;
818 static int nsets_in
=1,d
=1,nb_min
=4,resol
=10;
819 static real temp
=298.15,fit_start
=1,smooth_tail_start
=-1;
821 /* must correspond to enum avbar* declared at beginning of file */
822 static const char *avbar_opt
[avbarNR
+1] = {
823 NULL
, "none", "stddev", "error", "90", NULL
827 { "-time", FALSE
, etBOOL
, {&bHaveT
},
828 "Expect a time in the input" },
829 { "-b", FALSE
, etREAL
, {&tb
},
830 "First time to read from set" },
831 { "-e", FALSE
, etREAL
, {&te
},
832 "Last time to read from set" },
833 { "-n", FALSE
, etINT
, {&nsets_in
},
834 "Read # sets seperated by &" },
835 { "-d", FALSE
, etBOOL
, {&bDer
},
836 "Use the derivative" },
837 { "-dp", FALSE
, etINT
, {&d
},
838 "HIDDENThe derivative is the difference over # points" },
839 { "-bw", FALSE
, etREAL
, {&binwidth
},
840 "Binwidth for the distribution" },
841 { "-errbar", FALSE
, etENUM
, {avbar_opt
},
842 "Error bars for -av" },
843 { "-integrate",FALSE
,etBOOL
, {&bIntegrate
},
844 "Integrate data function(s) numerically using trapezium rule" },
845 { "-aver_start",FALSE
, etREAL
, {&aver_start
},
846 "Start averaging the integral from here" },
847 { "-xydy", FALSE
, etBOOL
, {&bXYdy
},
848 "Interpret second data set as error in the y values for integrating" },
849 { "-regression",FALSE
,etBOOL
,{&bRegression
},
850 "Perform a linear regression analysis on the data" },
851 { "-luzar", FALSE
, etBOOL
, {&bLuzar
},
852 "Do a Luzar and Chandler analysis on a correlation function and related as produced by g_hbond. When in addition the -xydy flag is given the second and fourth column will be interpreted as errors in c(t) and n(t)." },
853 { "-temp", FALSE
, etREAL
, {&temp
},
854 "Temperature for the Luzar hydrogen bonding kinetics analysis" },
855 { "-fitstart", FALSE
, etREAL
, {&fit_start
},
856 "Time (ps) from which to start fitting the correlation functions in order to obtain the forward and backward rate constants for HB breaking and formation" },
857 { "-smooth",FALSE
, etREAL
, {&smooth_tail_start
},
858 "If >= 0, the tail of the ACF will be smoothed by fitting it to an exponential function: y = A exp(-x/tau)" },
859 { "-nbmin", FALSE
, etINT
, {&nb_min
},
860 "HIDDENMinimum number of blocks for block averaging" },
861 { "-resol", FALSE
, etINT
, {&resol
},
862 "HIDDENResolution for the block averaging, block size increases with"
863 " a factor 2^(1/#)" },
864 { "-eeexpfit", FALSE
, etBOOL
, {&bEESEF
},
865 "HIDDENAlways use a single exponential fit for the error estimate" },
866 { "-eenlc", FALSE
, etBOOL
, {&bEENLC
},
867 "HIDDENAllow a negative long-time correlation" },
868 { "-eefitac", FALSE
, etBOOL
, {&bEeFitAc
},
869 "HIDDENAlso plot analytical block average using a autocorrelation fit" },
870 { "-filter", FALSE
, etREAL
, {&filtlen
},
871 "Print the high-frequency fluctuation after filtering with a cosine filter of length #" },
872 { "-power", FALSE
, etBOOL
, {&bPower
},
873 "Fit data to: b t^a" },
874 { "-subav", FALSE
, etBOOL
, {&bSubAv
},
875 "Subtract the average before autocorrelating" },
876 { "-oneacf", FALSE
, etBOOL
, {&bAverCorr
},
877 "Calculate one ACF over all sets" }
879 #define NPA asize(pa)
882 int n
,nlast
,s
,nset
,i
,j
=0;
883 real
**val
,*t
,dt
,tot
,error
;
884 double *av
,*sig
,cum1
,cum2
,cum3
,cum4
,db
;
885 const char *acfile
,*msdfile
,*ccfile
,*distfile
,*avfile
,*eefile
,*fitfile
;
889 { efXVG
, "-f", "graph", ffREAD
},
890 { efXVG
, "-ac", "autocorr", ffOPTWR
},
891 { efXVG
, "-msd", "msd", ffOPTWR
},
892 { efXVG
, "-cc", "coscont", ffOPTWR
},
893 { efXVG
, "-dist", "distr", ffOPTWR
},
894 { efXVG
, "-av", "average", ffOPTWR
},
895 { efXVG
, "-ee", "errest", ffOPTWR
},
896 { efLOG
, "-g", "fitlog", ffOPTWR
}
898 #define NFILE asize(fnm)
904 ppa
= add_acf_pargs(&npargs
,pa
);
906 CopyRight(stderr
,argv
[0]);
907 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
,
908 NFILE
,fnm
,npargs
,ppa
,asize(desc
),desc
,0,NULL
,&oenv
);
910 acfile
= opt2fn_null("-ac",NFILE
,fnm
);
911 msdfile
= opt2fn_null("-msd",NFILE
,fnm
);
912 ccfile
= opt2fn_null("-cc",NFILE
,fnm
);
913 distfile
= opt2fn_null("-dist",NFILE
,fnm
);
914 avfile
= opt2fn_null("-av",NFILE
,fnm
);
915 eefile
= opt2fn_null("-ee",NFILE
,fnm
);
916 if (opt2parg_bSet("-fitfn",npargs
,ppa
))
917 fitfile
= opt2fn("-g",NFILE
,fnm
);
919 fitfile
= opt2fn_null("-g",NFILE
,fnm
);
921 val
=read_xvg_time(opt2fn("-f",NFILE
,fnm
),bHaveT
,
922 opt2parg_bSet("-b",npargs
,ppa
),tb
,
923 opt2parg_bSet("-e",npargs
,ppa
),te
,
924 nsets_in
,&nset
,&n
,&dt
,&t
);
925 printf("Read %d sets of %d points, dt = %g\n\n",nset
,n
,dt
);
928 printf("Calculating the derivative as (f[i+%d]-f[i])/(%d*dt)\n\n",
931 for(s
=0; s
<nset
; s
++)
933 val
[s
][i
] = (val
[s
][i
+d
]-val
[s
][i
])/(d
*dt
);
937 printf("Calculating the integral using the trapezium rule\n");
940 sum
= evaluate_integral(n
,t
,val
[0],val
[1],aver_start
,&stddev
);
941 printf("Integral %10.3f +/- %10.5f\n",sum
,stddev
);
944 for(s
=0; s
<nset
; s
++) {
945 sum
= evaluate_integral(n
,t
,val
[s
],NULL
,aver_start
,&stddev
);
946 printf("Integral %d %10.5f +/- %10.5f\n",s
+1,sum
,stddev
);
951 out_fit
= ffopen(fitfile
,"w");
952 if (bXYdy
&& nset
>=2) {
953 do_fit(out_fit
,0,TRUE
,n
,t
,val
,npargs
,ppa
,oenv
);
955 for(s
=0; s
<nset
; s
++)
956 do_fit(out_fit
,s
,FALSE
,n
,t
,val
,npargs
,ppa
,oenv
);
961 printf(" std. dev. relative deviation of\n");
962 printf(" standard --------- cumulants from those of\n");
963 printf("set average deviation sqrt(n-1) a Gaussian distribition\n");
964 printf(" cum. 3 cum. 4\n");
967 for(s
=0; (s
<nset
); s
++) {
975 for(i
=0; (i
<n
); i
++) {
987 error
= sqrt(cum2
/(n
-1));
990 printf("SS%d %13.6e %12.6e %12.6e %6.3f %6.3f\n",
991 s
+1,av
[s
],sig
[s
],error
,
992 sig
[s
] ? cum3
/(sig
[s
]*sig
[s
]*sig
[s
]*sqrt(8/M_PI
)) : 0,
993 sig
[s
] ? cum4
/(sig
[s
]*sig
[s
]*sig
[s
]*sig
[s
]*3)-1 : 0);
998 filter(filtlen
,n
,nset
,val
,dt
,oenv
);
1001 out
=xvgropen(msdfile
,"Mean square displacement",
1002 "time","MSD (nm\\S2\\N)",oenv
);
1003 nlast
= (int)(n
*frac
);
1004 for(s
=0; s
<nset
; s
++) {
1005 for(j
=0; j
<=nlast
; j
++) {
1007 fprintf(stderr
,"\r%d",j
);
1009 for(i
=0; i
<n
-j
; i
++)
1010 tot
+= sqr(val
[s
][i
]-val
[s
][i
+j
]);
1012 fprintf(out
," %g %8g\n",dt
*j
,tot
);
1018 fprintf(stderr
,"\r%d, time=%g\n",j
-1,(j
-1)*dt
);
1021 plot_coscont(ccfile
,n
,nset
,val
,oenv
);
1024 histogram(distfile
,binwidth
,n
,nset
,val
,oenv
);
1026 average(avfile
,nenum(avbar_opt
),n
,nset
,val
,t
);
1028 estimate_error(eefile
,nb_min
,resol
,n
,nset
,av
,sig
,val
,dt
,
1029 bEeFitAc
,bEESEF
,bEENLC
,oenv
);
1031 power_fit(n
,nset
,val
,t
);
1034 for(s
=0; s
<nset
; s
++)
1037 do_autocorr(acfile
,oenv
,"Autocorrelation",n
,nset
,val
,dt
,
1038 eacNormal
,bAverCorr
);
1041 regression_analysis(n
,bXYdy
,t
,val
);
1044 luzar_correl(n
,t
,nset
,val
,temp
,bXYdy
,fit_start
,smooth_tail_start
,oenv
);
1046 view_all(oenv
,NFILE
, fnm
);