File: | gromacs/gmxana/autocorr.c |
Location: | line 905, column 13 |
Description: | Value stored to 'sum' is never read |
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 | |
59 | typedef struct { |
60 | unsigned long mode; |
61 | int nrestart, nout, P, fitfn; |
62 | gmx_bool bFour, bNormalize; |
63 | real tbeginfit, tendfit; |
64 | } t_acf; |
65 | |
66 | static gmx_bool bACFinit = FALSE0; |
67 | static t_acf acf; |
68 | |
69 | enum { |
70 | enNorm, enCos, enSin |
71 | }; |
72 | |
73 | int 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 | |
89 | static 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 | |
156 | static 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 | |
262 | void 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 | |
302 | void 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 | |
323 | void 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 | |
339 | void 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 | |
352 | real 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 | |
389 | real 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 | |
445 | void 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 | |
636 | real 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]; |
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 | |
776 | void 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); |
Value stored to 'sum' is never read | |
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 | |
972 | static const char *Leg[] = { NULL((void*)0), "0", "1", "2", "3", NULL((void*)0) }; |
973 | |
974 | t_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 | |
1024 | void 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 | |
1058 | int 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 | |
1068 | int 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 | |
1078 | void 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 | } |