Bug Summary

File:gromacs/gmxana/autocorr.c
Location:line 664, column 5
Description:Value stored to 'ct_estimate' is never read

Annotated Source Code

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