?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / tools / autocorr.c
blob8c35084c75148741b5c724e9fec3ab8d81301eab
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_autocorr_c = "$Id$";
30
31 #include <stdio.h>
32 #include <math.h>
33 #include "macros.h"
34 #include "typedefs.h"
35 #include "physics.h"
36 #include "smalloc.h"
37 #include "xvgr.h"
38 #include "futil.h"
39 #include "gstat.h"
40 #include "names.h"
41 #include "fatal.h"
42 #include "vec.h"
43 #include "string2.h"
44
45 #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
46 #define MODE(x) ((mode & (x)) == (x))
47
48 typedef struct {
49 unsigned long mode;
50 int nrestart,nout,P,nfitparm;
51 bool bFour,bNormalize;
52 real tbeginfit,tendfit;
53 } t_acf;
54
55 static bool bACFinit = FALSE;
56 static t_acf acf;
57
58 typedef real fftreal;
59
60 void four1(fftreal data[],int nn,int isign)
61 {
62 int n,mmax,m,j,istep,i;
63 double wtemp,wr,wpr,wpi,wi,theta;
64 fftreal tempr,tempi;
65
66 n=nn << 1;
67 j=1;
68 for (i=1;i<n;i+=2) {
69 if (j > i) {
70 SWAP(data[j],data[i]);
71 SWAP(data[j+1],data[i+1]);
72 }
73 m=n >> 1;
74 while (m >= 2 && j > m) {
75 j -= m;
76 m >>= 1;
77 }
78 j += m;
79 }
80 mmax=2;
81 while (n > mmax) {
82 istep=2*mmax;
83 theta=6.28318530717959/(isign*mmax);
84 wtemp=sin(0.5*theta);
85 wpr = -2.0*wtemp*wtemp;
86 wpi=sin(theta);
87 wr=1.0;
88 wi=0.0;
89 for (m=1;m<mmax;m+=2) {
90 for (i=m;i<=n;i+=istep) {
91 j=i+mmax;
92 tempr=wr*data[j]-wi*data[j+1];
93 tempi=wr*data[j+1]+wi*data[j];
94 data[j]=data[i]-tempr;
95 data[j+1]=data[i+1]-tempi;
96 data[i] += tempr;
97 data[i+1] += tempi;
98 }
99 wr=(wtemp=wr)*wpr-wi*wpi+wr;
100 wi=wi*wpr+wtemp*wpi+wi;
101 }
102 mmax=istep;
103 }
104 }
105
106 #undef SWAP
107
108 void realft(fftreal data[],int n,int isign)
109 {
110 int i,i1,i2,i3,i4,n2p3;
111 fftreal c1=0.5,c2,h1r,h1i,h2r,h2i;
112 double wr,wi,wpr,wpi,wtemp,theta;
113
114 theta=3.141592653589793/(double) n;
115 if (isign == 1) {
116 c2 = -0.5;
117 four1(data,n,1);
118 } else {
119 c2=0.5;
120 theta = -theta;
121 }
122 wtemp=sin(0.5*theta);
123 wpr = -2.0*wtemp*wtemp;
124 wpi=sin(theta);
125 wr=1.0+wpr;
126 wi=wpi;
127 n2p3=2*n+3;
128 for (i=2;i<=n/2;i++) {
129 i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
130 h1r=c1*(data[i1]+data[i3]);
131 h1i=c1*(data[i2]-data[i4]);
132 h2r = -c2*(data[i2]+data[i4]);
133 h2i=c2*(data[i1]-data[i3]);
134 data[i1]=h1r+wr*h2r-wi*h2i;
135 data[i2]=h1i+wr*h2i+wi*h2r;
136 data[i3]=h1r-wr*h2r+wi*h2i;
137 data[i4] = -h1i+wr*h2i+wi*h2r;
138 wr=(wtemp=wr)*wpr-wi*wpi+wr;
139 wi=wi*wpr+wtemp*wpi+wi;
140 }
141 if (isign == 1) {
142 data[1] = (h1r=data[1])+data[2];
143 data[2] = h1r-data[2];
144 } else {
145 data[1]=c1*((h1r=data[1])+data[2]);
146 data[2]=c1*(h1r-data[2]);
147 four1(data,n,-1);
148 }
149 }
150
151 void twofft(fftreal data1[],fftreal data2[],fftreal fft1[],fftreal fft2[],int n)
152 {
153 int nn3,nn2,jj,j;
154 fftreal rep,rem,aip,aim;
155
156 nn3=1+(nn2=2+n+n);
157 for (j=1,jj=2;j<=n;j++,jj+=2) {
158 fft1[jj-1]=data1[j];
159 fft1[jj]=data2[j];
160 }
161 four1(fft1,n,1);
162 fft2[1]=fft1[2];
163 fft1[2]=fft2[2]=0.0;
164 for (j=3;j<=n+1;j+=2) {
165 rep=0.5*(fft1[j]+fft1[nn2-j]);
166 rem=0.5*(fft1[j]-fft1[nn2-j]);
167 aip=0.5*(fft1[j+1]+fft1[nn3-j]);
168 aim=0.5*(fft1[j+1]-fft1[nn3-j]);
169 fft1[j]=rep;
170 fft1[j+1]=aim;
171 fft1[nn2-j]=rep;
172 fft1[nn3-j] = -aim;
173 fft2[j]=aip;
174 fft2[j+1] = -rem;
175 fft2[nn2-j]=aip;
176 fft2[nn3-j]=rem;
177 }
178 }
179
180 enum { enNorm, enCos, enSin };
181
182 void correl(fftreal data1[],fftreal data2[],int n,fftreal ans[])
183 {
184 int no2,i;
185 fftreal dum,*fft;
186
187 snew(fft,2*n+1);
188 twofft(data1,data2,fft,ans,n);
189 no2=n/2;
190 for (i=2;i<=n+2;i+=2) {
191 dum = ans[i-1];
192 ans[i-1] = (fft[i-1]*dum+fft[i]*ans[i])/no2;
193 ans[i] = (fft[i]*dum-fft[i-1]*ans[i])/no2;
194 }
195 ans[2]=ans[n+1];
196 realft(ans,no2,-1);
197 sfree(fft);
198 }
199
200 static void low_do_four_core(int nfour,int nframes,real c1[],fftreal cfour[],
201 int nCos,bool bPadding)
202 {
203 int i=0;
204 fftreal aver,*ans;
205
206 aver = 0.0;
207 switch (nCos) {
208 case enNorm:
209 for(i=0; (i<nframes); i++) {
210 aver+=c1[i];
211 cfour[i]=c1[i];
212 }
213 break;
214 case enCos:
215 for(i=0; (i<nframes); i++)
216 cfour[i]=cos(c1[i]);
217 break;
218 case enSin:
219 for(i=0; (i<nframes); i++)
220 cfour[i]=sin(c1[i]);
221 break;
222 default:
223 fatal_error(0,"nCos = %d, %s %d",nCos,__FILE__,__LINE__);
224 }
225 for( ; (i<nfour); i++)
226 cfour[i]= 0.0;
227
228 if (bPadding) {
229 aver /= nframes;
230 /* fprintf(stderr,"aver = %g\n",aver); */
231 for(i=0; (i<nframes); i++)
232 cfour[i] -= aver;
233 }
234
235 snew(ans,2*nfour);
236 correl(cfour-1,cfour-1,nfour,ans-1);
237
238 if (bPadding)
239 for (i=0; (i<nfour); i++)
240 cfour[i] = ans[i]+sqr(aver);
241 else
242 for (i=0; (i<nfour); i++)
243 cfour[i] = ans[i];
244
245 sfree(ans);
246 }
247
248 static void do_ac_core(int nframes,int nout,
249 real corr[],real c1[],int nrestart,
250 unsigned long mode)
251 {
252 int j,k,j3,jk3,m,n;
253 fftreal ccc,c0,cth;
254 rvec xj,xk,rr;
255
256 if (nrestart < 1) {
257 fprintf(stderr,"WARNING: setting number of restarts to 1\n");
258 nrestart = 1;
259 }
260 if (debug)
261 fprintf(debug,
262 "Starting do_ac_core: nframes=%d, nout=%d, nrestart=%d,mode=%lu\n",
263 nframes,nout,nrestart,mode);
264
265 for(j=0; (j<nout); j++)
266 corr[j]=0;
267
268 /* Loop over starting points. */
269 for(j=0; (j<nframes); j+=nrestart) {
270 j3 = DIM*j;
271
272 /* Loop over the correlation length for this starting point */
273 for(k=0; (k<nout) && (j+k < nframes); k++) {
274 jk3 = DIM*(j+k);
275
276 /* Switch over possible ACF types.
277 * It might be more efficient to put the loops inside the switch,
278 * but this is more clear, and save development time!
279 */
280 if (MODE(eacNormal)) {
281 corr[k] += c1[j]*c1[j+k];
282 }
283 else if (MODE(eacCos)) {
284 /* Compute the cos (phi(t)-phi(t+dt)) */
285 corr[k] += cos(c1[j]-c1[j+k]);
286 }
287 else if (MODE(eacP1) || MODE(eacP2) || MODE(eacP3)) {
288 for(m=0; (m<DIM); m++) {
289 xj[m] = c1[j3+m];
290 xk[m] = c1[jk3+m];
291 }
292 cth=cos_angle(xj,xk);
293
294 if (cth-1.0 > 1.0e-15) {
295 fprintf(stderr,"j: %d, k: %d, xj:(%g,%g,%g), xk:(%g,%g,%g)\n",
296 j,k,xj[XX],xj[YY],xj[ZZ],xk[XX],xk[YY],xk[ZZ]);
297 }
298
299 corr[k] += LegendreP(cth,mode); /* 1.5*cth*cth-0.5; */
300 }
301 else if (MODE(eacRcross)) {
302 for(m=0; (m<DIM); m++) {
303 xj[m] = c1[j3+m];
304 xk[m] = c1[jk3+m];
305 }
306 oprod(xj,xk,rr);
307
308 corr[k] += iprod(rr,rr);
309 }
310 else if (MODE(eacVector)) {
311 for(m=0; (m<DIM); m++) {
312 xj[m] = c1[j3+m];
313 xk[m] = c1[jk3+m];
314 }
315 ccc = iprod(xj,xk);
316
317 corr[k] += ccc;
318 }
319 else
320 fatal_error(0,"\nInvalid mode (%d) in do_ac_core",mode);
321 }
322 }
323 /* Correct for the number of points and copy results to the data array */
324 for(j=0; (j<nout); j++) {
325 n = (nframes-j+(nrestart-1))/nrestart;
326 c1[j] = corr[j]/n;
327 }
328 }
329
330 void normalize_acf(int nout,real corr[])
331 {
332 int j;
333 real c0;
334
335 if (debug) {
336 fprintf(debug,"Before normalization\n");
337 for(j=0; (j<nout); j++)
338 fprintf(debug,"%5d %10f\n",j,corr[j]);
339 }
340
341 /* Normalisation makes that c[0] = 1.0 and that other points are scaled
342 * accordingly.
343 */
344 if (corr[0] == 0.0)
345 c0 = 1.0;
346 else
347 c0 = 1.0/corr[0];
348 for(j=0; (j<nout); j++)
349 corr[j] *= c0;
350
351 if (debug) {
352 fprintf(debug,"After normalization\n");
353 for(j=0; (j<nout); j++)
354 fprintf(debug,"%5d %10f\n",j,corr[j]);
355 }
356 }
357
358 void average_acf(int n,int nitem,real **c1)
359 {
360 real c0;
361 int i,j;
362
363 fprintf(stderr,"Averaging correlation functions\n");
364
365 for(j=0; (j<n); j++) {
366 c0 = 0;
367 for(i=0; (i<nitem); i++)
368 c0+=c1[i][j];
369 c1[0][j] = c0/nitem;
370 }
371 }
372
373 void norm_and_scale_vectors(int nframes,real c1[],real scale)
374 {
375 int j,m;
376 real *rij;
377
378 for(j=0; (j<nframes); j++) {
379 rij = &(c1[j*DIM]);
380 unitv(rij,rij);
381 for(m=0; (m<DIM); m++)
382 rij[m]*=scale;
383 }
384 }
385
386 void dump_tmp(char *s,int n,real c[])
387 {
388 FILE *fp;
389 int i;
390
391 fp=ffopen(s,"w");
392 for(i=0; (i<n); i++)
393 fprintf(fp,"%10d %10g\n",i,c[i]);
394 ffclose(fp);
395 }
396
397 real print_and_integrate(FILE *fp,int n,real dt,real c[])
398 {
399 real c0,sum;
400 int j;
401
402 /* Use trapezoidal rule for calculating integral */
403 sum = 0.0;
404 for(j=0; (j<n); j++) {
405 c0 = c[j];
406 if (fp)
407 fprintf(fp,"%10.3f %10.5f\n",j*dt,c0);
408 if (j > 0)
409 sum+=dt*(c0+c[j-1]);
410 }
411 fprintf(fp,"&\n");
412 return sum*0.5;
413 }
414
415 void do_four_core(unsigned long mode,int nfour,int nf2,int nframes,
416 real c1[],real csum[],real ctmp[])
417 {
418 fftreal *cfour;
419 char buf[32];
420 real fac;
421 int i,j,m,m1;
422
423 snew(cfour,nfour);
424
425 if (MODE(eacNormal)) {
426 /********************************************
427 * N O R M A L
428 ********************************************/
429 low_do_four_core(nfour,nf2,c1,csum,enNorm,FALSE);
430 }
431 else if (MODE(eacCos)) {
432 /***************************************************
433 * C O S I N E
434 ***************************************************/
435 /* Copy the data to temp array. Since we need it twice
436 * we can't overwrite original.
437 */
438 for(j=0; (j<nf2); j++)
439 ctmp[j]=c1[j];
440
441 /* Cosine term of AC function */
442 low_do_four_core(nfour,nf2,ctmp,cfour,enCos,FALSE);
443 for(j=0; (j<nf2); j++)
444 c1[j] = cfour[j];
445
446 /* Sine term of AC function */
447 low_do_four_core(nfour,nf2,ctmp,cfour,enSin,FALSE);
448 for(j=0; (j<nf2); j++) {
449 c1[j] += cfour[j];
450 csum[j] = c1[j];
451 }
452 }
453 else if (MODE(eacP2)) {
454 /***************************************************
455 * Legendre polynomials
456 ***************************************************/
457 /* First normalize the vectors */
458 norm_and_scale_vectors(nframes,c1,1.0);
459
460 /* For P2 thingies we have to do six FFT based correls
461 * First for XX^2, then for YY^2, then for ZZ^2
462 * Then we have to do XY, YZ and XZ (counting these twice)
463 * After that we sum them and normalise
464 * P2(x) = (3 * cos^2 (x) - 1)/2
465 * for unit vectors u and v we compute the cosine as the inner product
466 * cos(u,v) = uX vX + uY vY + uZ vZ
467 *
468 * oo
469 * /
470 * C(t) = | (3 cos^2(u(t'),u(t'+t)) - 1)/2 dt'
471 * /
472 * 0
473 *
474 * For ACF we need:
475 * P2(u(0),u(t)) = [3 * (uX(0) uX(t) +
476 * uY(0) uY(t) +
477 * uZ(0) uZ(t))^2 - 1]/2
478 * = [3 * ((uX(0) uX(t))^2 +
479 * (uY(0) uY(t))^2 +
480 * (uZ(0) uZ(t))^2 +
481 * 2(uX(0) uY(0) uX(t) uY(t)) +
482 * 2(uX(0) uZ(0) uX(t) uZ(t)) +
483 * 2(uY(0) uZ(0) uY(t) uZ(t))) - 1]/2
484 *
485 * = [(3/2) * (<uX^2> + <uY^2> + <uZ^2> +
486 * 2<uXuY> + 2<uXuZ> + 2<uYuZ>) - 0.5]
487 *
488 */
489
490 /* Because of normalization the number of -0.5 to subtract
491 * depends on the number of data points!
492 */
493 for(j=0; (j<nf2); j++)
494 csum[j] = -0.5*(nf2-j);
495
496 /***** DIAGONAL ELEMENTS ************/
497 for(m=0; (m<DIM); m++) {
498 /* Copy the vector data in a linear array */
499 for(j=0; (j<nf2); j++)
500 ctmp[j] = sqr(c1[DIM*j+m]);
501 if (debug) {
502 sprintf(buf,"c1diag%d.xvg",m);
503 dump_tmp(buf,nf2,ctmp);
504 }
505
506 low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
507
508 if (debug) {
509 sprintf(buf,"c1dfout%d.xvg",m);
510 dump_tmp(buf,nf2,cfour);
511 }
512 fac = 1.5;
513 for(j=0; (j<nf2); j++)
514 csum[j] += fac*(cfour[j]);
515 }
516 /******* OFF-DIAGONAL ELEMENTS **********/
517 for(m=0; (m<DIM); m++) {
518 /* Copy the vector data in a linear array */
519 m1=(m+1) % DIM;
520 for(j=0; (j<nf2); j++)
521 ctmp[j]=c1[DIM*j+m]*c1[DIM*j+m1];
522
523 if (debug) {
524 sprintf(buf,"c1off%d.xvg",m);
525 dump_tmp(buf,nf2,ctmp);
526 }
527 low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
528 if (debug) {
529 sprintf(buf,"c1ofout%d.xvg",m);
530 dump_tmp(buf,nf2,cfour);
531 }
532 fac = 3.0;
533 for(j=0; (j<nf2); j++) {
534 csum[j] += fac*cfour[j];
535 }
536 }
537 }
538 else if (MODE(eacP1) || MODE(eacVector)) {
539 /***************************************************
540 * V E C T O R & P1
541 ***************************************************/
542 if (MODE(eacP1)) {
543 /* First normalize the vectors */
544 norm_and_scale_vectors(nframes,c1,1.0);
545 }
546
547 /* For vector thingies we have to do three FFT based correls
548 * First for XX, then for YY, then for ZZ
549 * After that we sum them and normalise
550 */
551 for(j=0; (j<nf2); j++) {
552 csum[j]=0.0;
553 }
554 for(m=0; (m<DIM); m++) {
555 /* Copy the vector data in a linear array */
556 for(j=0; (j<nf2); j++)
557 ctmp[j]=c1[DIM*j+m];
558 low_do_four_core(nfour,nf2,ctmp,cfour,enNorm,FALSE);
559 for(j=0; (j<nf2); j++)
560 csum[j] += cfour[j];
561 }
562 }
563 else
564 fatal_error(0,"\nUnknown mode in do_autocorr (%d)",mode);
565
566 sfree(cfour);
567 for(j=0; (j<nf2); j++)
568 c1[j] = csum[j]/(real)(nframes-j);
569 }
570
571 void fit_acf(int ncorr,int nfitparm,bool bVerbose,
572 real tbeginfit,real tendfit,real dt,real c1[])
573 {
574 real fitparm[3];
575 real tStart,tail_corr,sum,sumtot,*sig;
576 int j,nf_int;
577
578 fprintf(stderr,"CORR:\n");
579
580 nf_int = min(ncorr,(int)(tendfit/dt));
581 sum = print_and_integrate(debug,ncorr,dt,c1);
582
583 fprintf(stderr,"CORR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
584 0.0,dt*nf_int,sum);
585 fprintf(stderr,"CORR: Relaxation times are computed as fit to an exponential:\n");
586 if (nfitparm == 1)
587 fprintf(stderr,"CORR: Exp[-t/tau_slope]\n");
588 else if (nfitparm == 2)
589 fprintf(stderr,"CORR: A Exp[-t/tau_slope]\n");
590 else
591 fatal_error(0,"nparm not set to 1 or 2, %s %d",__FILE__,__LINE__);
592 fprintf(stderr,"CORR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n",tbeginfit,min(ncorr*dt,tendfit));
593
594 tStart = 0;
595 fprintf(stderr,"CORR:%12s%12s%12s%12s%12s\n",
596 "Integral to","Value","Tail Value","Sum (ps)","Tau (ps)");
597 for(j=0; ((j<5) && (tStart < tendfit)); j++) {
598 snew(sig,ncorr);
599 fitparm[0]=fitparm[1]=fitparm[2] = 1.0;
600 nf_int = min(ncorr,(int)((tStart+1e-4)/dt));
601 sum = print_and_integrate(debug,nf_int,dt,c1);
602 tail_corr = do_lmfit(ncorr,c1,sig,dt,tStart,tendfit,
603 bVerbose,nfitparm,NULL,fitparm,NULL);
604 sumtot = sum+tail_corr;
605 fprintf(stderr,"CORR:%12.5e%12.5e%12.5e%12.5e%12.5e\n",
606 tStart,sum,tail_corr,sumtot,fitparm[0]);
607
608 tStart += tbeginfit;
609 sfree(sig);
610 }
611 }
612
613 void low_do_autocorr(char *fn,char *title,
614 int nframes,int nitem,int nout,real **c1,
615 real dt,unsigned long mode,int nrestart,
616 bool bAver,bool bFour,bool bNormalize,
617 bool bVerbose,real tbeginfit,real tendfit,
618 int nfitparm)
619 {
620 FILE *fp;
621 int i,k,nfour;
622 fftreal *csum;
623 real *ctmp;
624 real c0,sum;
625
626 /* Check flags and parameters */
627 /* nout = get_acfnout();*/
628 if (nout == -1)
629 nout = acf.nout = (nframes+1)/2;
630 else if (nout > nframes)
631 nout=nframes;
632
633 if (MODE(eacCos) && MODE(eacVector))
634 fatal_error(0,"Incompatible options bCos && bVector (%s, %d)",
635 __FILE__,__LINE__);
636 if ((MODE(eacP3) || MODE(eacRcross)) && bFour) {
637 fprintf(stderr,"Can't combine mode %lu with FFT, turning off FFT\n",mode);
638 bFour = FALSE;
639 }
640 if (MODE(eacNormal) && MODE(eacVector))
641 fatal_error(0,"Incompatible mode bits: normal and vector (or Legendre)");
642
643 /* Print flags and parameters */
644 fprintf(stderr,"Will calculate %s of %d thingies for %d frames\n",
645 title,nitem,nframes);
646 fprintf(stderr,"bAver = %s, bFour = %s bNormalize= %s\n",
647 bool_names[bAver],bool_names[bFour],bool_names[bNormalize]);
648 fprintf(stderr,"mode = %lu, dt = %g, nrestart = %d\n",mode,dt,nrestart);
649
650 if (bFour) {
651 c0 = log((double)nframes)/log(2.0);
652 k = c0;
653 if (k < c0)
654 k++;
655 k++;
656 nfour = 1<<k;
657 if (debug)
658 fprintf(debug,"Using FFT to calculate %s, #points for FFT = %d\n",
659 title,nfour);
660
661 /* Allocate temp arrays */
662 snew(csum,nfour);
663 snew(ctmp,nfour);
664 }
665 else {
666 nfour = 0; /* To keep the compiler happy */
667 snew(csum,nframes);
668 snew(ctmp,nframes);
669 }
670
671 /* Loop over items (e.g. molecules or dihedrals)
672 * In this loop the actual correlation functions are computed, but without
673 * normalizing them.
674 */
675 for(i=0; (i<nitem); i++) {
676 fprintf(stderr,"\rThingie %d",i);
677
678 if (bFour)
679 do_four_core(mode,nfour,nframes,nframes,c1[i],csum,ctmp);
680 else
681 do_ac_core(nframes,nout,ctmp,c1[i],nrestart,mode);
682 }
683 fprintf(stderr,"\n");
684 sfree(ctmp);
685 sfree(csum);
686
687 fp=xvgropen(fn,title,"Time (ps)","C(t)");
688 if (bAver) {
689 average_acf(nframes,nitem,c1);
690
691 if (bNormalize)
692 normalize_acf(nout,c1[0]);
693
694 if (tbeginfit < tendfit) {
695 fit_acf(nout,nfitparm,bVerbose,tbeginfit,tendfit,dt,c1[0]);
696 (void)print_and_integrate(fp,nout,dt,c1[0]);
697 } else {
698 sum = print_and_integrate(fp,nout,dt,c1[0]);
699 fprintf(stderr,"Correlation time (integral over corrfn): %g (ps)\n",sum);
700 }
701 }
702 else {
703 /* Not averaging. Normalize individual ACFs */
704 for(i=0; (i<nitem); i++) {
705 if (bNormalize)
706 normalize_acf(nout,c1[i]);
707 if (tbeginfit < tendfit) {
708 fit_acf(nout,nfitparm,bVerbose,tbeginfit,tendfit,dt,c1[i]);
709 (void)print_and_integrate(fp,nout,dt,c1[0]);
710 } else {
711 sum = print_and_integrate(fp,nout,dt,c1[i]);
712 fprintf(stderr,"CORRelation time (integral over corrfn %d): %g (ps)\n",
713 i,sum);
714 }
715 }
716 }
717 ffclose(fp);
718 }
719
720 static char *Leg[] = { NULL, "0", "1", "2", "3", NULL };
721 static char *Nparm[] = { NULL, "1", "2", NULL };
722
723 t_pargs *add_acf_pargs(int *npargs,t_pargs *pa)
724 {
725 t_pargs acfpa[] = {
726 { "-acflen", FALSE, etINT, {&acf.nout},
727 "Length of the ACF, default is half the number of frames" },
728 { "-normalize",FALSE, etBOOL, {&acf.bNormalize},
729 "Normalize ACF" },
730 { "-fft", FALSE, etBOOL, {&acf.bFour},
731 "HIDDENUse fast fourier transform for correlation function" },
732 { "-nrestart", FALSE, etINT, {&acf.nrestart},
733 "HIDDENNumber of frames between time origins for ACF when no FFT is used" },
734 { "-P", FALSE, etENUM, {Leg},
735 "Order of Legendre polynomial for ACF (0 indicates none)" },
736 { "-nparm", FALSE, etENUM, {Nparm},
737 "Number of parameters in exponential fit" },
738 { "-beginfit", FALSE, etREAL, {&acf.tbeginfit},
739 "Time where to begin the exponential fit of the correlation function" },
740 { "-endfit", FALSE, etREAL, {&acf.tendfit},
741 "Time where to end the exponential fit of the correlation function" },
742 };
743 #define NPA asize(acfpa)
744 t_pargs *ppa;
745 int i;
746
747 snew(ppa,*npargs+NPA);
748 for(i=0; (i<*npargs); i++)
749 ppa[i] = pa[i];
750 for(i=0; (i<NPA); i++)
751 ppa[*npargs+i] = acfpa[i];
752 (*npargs) += NPA;
753
754 acf.mode = 0;
755 acf.nrestart = 1;
756 acf.nout = -1;
757 acf.P = 0;
758 acf.nfitparm = 1;
759 acf.bFour = TRUE;
760 acf.bNormalize = TRUE;
761 acf.tbeginfit = 0.0;
762 acf.tendfit = 0.0;
763
764 bACFinit = TRUE;
765
766 return ppa;
767 }
768
769 void do_autocorr(char *fn,char *title,int nframes,int nitem,real **c1,
770 real dt,unsigned long mode,bool bAver)
771 {
772 if (!bACFinit) {
773 fprintf(stderr,"ACF data structures have not been initialised. Call add_acf_pargs\n");
774 }
775
776 /* Handle enumerated types */
777 sscanf(Leg[0],"%d",&acf.P);
778 sscanf(Nparm[0],"%d",&acf.nfitparm);
779
780 switch (acf.P) {
781 case 1:
782 mode = mode | eacP1;
783 break;
784 case 2:
785 mode = mode | eacP2;
786 break;
787 case 3:
788 mode = mode | eacP3;
789 break;
790 default:
791 break;
792 }
793
794 low_do_autocorr(fn,title,nframes,nitem,acf.nout,c1,dt,mode,
795 acf.nrestart,bAver,acf.bFour,acf.bNormalize,
796 bDebugMode(),acf.tbeginfit,acf.tendfit,
797 acf.nfitparm);
798 }
799
800 int get_acfnout(void)
801 {
802 if (!bACFinit)
803 fatal_error(0,"ACF data not initialized yet");
804
805 return acf.nout;
806 }
got merge a long time ago