Bug Summary

File:gromacs/gmxana/anadih.c
Location:line 513, column 17
Description:Value stored to 'b' 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#include "physics.h"
45#include "gromacs/utility/smalloc.h"
46#include "macros.h"
47#include "txtdump.h"
48#include "bondf.h"
49#include "gromacs/fileio/xvgr.h"
50#include "typedefs.h"
51#include "gromacs/math/vec.h"
52#include "gstat.h"
53#include "gromacs/fileio/confio.h"
54#include "gromacs/fileio/trxio.h"
55
56#include "gromacs/utility/fatalerror.h"
57
58void print_one(const output_env_t oenv, const char *base, const char *name,
59 const char *title, const char *ylabel, int nf, real time[],
60 real data[])
61{
62 FILE *fp;
63 char buf[256], t2[256];
64 int k;
65
66 sprintf(buf, "%s%s.xvg", base, name);
67 fprintf(stderrstderr, "\rPrinting %s ", buf);
68 sprintf(t2, "%s %s", title, name);
69 fp = xvgropen(buf, t2, "Time (ps)", ylabel, oenv);
70 for (k = 0; (k < nf); k++)
71 {
72 fprintf(fp, "%10g %10g\n", time[k], data[k]);
73 }
74 gmx_ffclose(fp);
75}
76
77static int calc_RBbin(real phi, int gmx_unused__attribute__ ((unused)) multiplicity, real gmx_unused__attribute__ ((unused)) core_frac)
78{
79 /* multiplicity and core_frac NOT used,
80 * just given to enable use of pt-to-fn in caller low_ana_dih_trans*/
81 static const real r30 = M_PI3.14159265358979323846/6.0;
82 static const real r90 = M_PI3.14159265358979323846/2.0;
83 static const real r150 = M_PI3.14159265358979323846*5.0/6.0;
84
85 if ((phi < r30) && (phi > -r30))
86 {
87 return 1;
88 }
89 else if ((phi > -r150) && (phi < -r90))
90 {
91 return 2;
92 }
93 else if ((phi < r150) && (phi > r90))
94 {
95 return 3;
96 }
97 return 0;
98}
99
100static int calc_Nbin(real phi, int multiplicity, real core_frac)
101{
102 static const real r360 = 360*DEG2RAD(3.14159265358979323846/180.0);
103 real rot_width, core_width, core_offset, low, hi;
104 int bin;
105 /* with multiplicity 3 and core_frac 0.5
106 * 0<g(-)<120, 120<t<240, 240<g(+)<360
107 * 0< bin0 < 30, 30<bin1<90, 90<bin0<150, 150<bin2<210, 210<bin0<270, 270<bin3<330, 330<bin0<360
108 * so with multiplicity 3, bin1 is core g(-), bin2 is core t, bin3 is
109 core g(+), bin0 is between rotamers */
110 if (phi < 0)
111 {
112 phi += r360;
113 }
114
115 rot_width = 360/multiplicity;
116 core_width = core_frac * rot_width;
117 core_offset = (rot_width - core_width)/2.0;
118 for (bin = 1; bin <= multiplicity; bin++)
119 {
120 low = ((bin - 1) * rot_width ) + core_offset;
121 hi = ((bin - 1) * rot_width ) + core_offset + core_width;
122 low *= DEG2RAD(3.14159265358979323846/180.0);
123 hi *= DEG2RAD(3.14159265358979323846/180.0);
124 if ((phi > low) && (phi < hi))
125 {
126 return bin;
127 }
128 }
129 return 0;
130}
131
132void ana_dih_trans(const char *fn_trans, const char *fn_histo,
133 real **dih, int nframes, int nangles,
134 const char *grpname, real *time, gmx_bool bRb,
135 const output_env_t oenv)
136{
137 /* just a wrapper; declare extra args, then chuck away at end. */
138 int maxchi = 0;
139 t_dlist *dlist;
140 int *multiplicity;
141 int nlist = nangles;
142 int k;
143
144 snew(dlist, nlist)(dlist) = save_calloc("dlist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 144, (nlist), sizeof(*(dlist)));
145 snew(multiplicity, nangles)(multiplicity) = save_calloc("multiplicity", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 145, (nangles), sizeof(*(multiplicity)));
146 for (k = 0; (k < nangles); k++)
147 {
148 multiplicity[k] = 3;
149 }
150
151 low_ana_dih_trans(TRUE1, fn_trans, TRUE1, fn_histo, maxchi,
152 dih, nlist, dlist, nframes,
153 nangles, grpname, multiplicity, time, bRb, 0.5, oenv);
154 sfree(dlist)save_free("dlist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 154, (dlist));
155 sfree(multiplicity)save_free("multiplicity", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 155, (multiplicity));
156
157}
158
159void low_ana_dih_trans(gmx_bool bTrans, const char *fn_trans,
160 gmx_bool bHisto, const char *fn_histo, int maxchi,
161 real **dih, int nlist, t_dlist dlist[], int nframes,
162 int nangles, const char *grpname, int multiplicity[],
163 real *time, gmx_bool bRb, real core_frac,
164 const output_env_t oenv)
165{
166 FILE *fp;
167 int *tr_f, *tr_h;
168 char title[256];
169 int i, j, k, Dih, ntrans;
170 int cur_bin, new_bin;
171 real ttime, tt;
172 real *rot_occ[NROT4];
173 int (*calc_bin)(real, int, real);
174 real dt;
175
176 if (1 <= nframes)
177 {
178 return;
179 }
180 /* Assumes the frames are equally spaced in time */
181 dt = (time[nframes-1]-time[0])/(nframes-1);
182
183 /* Analysis of dihedral transitions */
184 fprintf(stderrstderr, "Now calculating transitions...\n");
185
186 if (bRb)
187 {
188 calc_bin = calc_RBbin;
189 }
190 else
191 {
192 calc_bin = calc_Nbin;
193 }
194
195 for (k = 0; k < NROT4; k++)
196 {
197 snew(rot_occ[k], nangles)(rot_occ[k]) = save_calloc("rot_occ[k]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 197, (nangles), sizeof(*(rot_occ[k])));
198 for (i = 0; (i < nangles); i++)
199 {
200 rot_occ[k][i] = 0;
201 }
202 }
203 snew(tr_h, nangles)(tr_h) = save_calloc("tr_h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 203, (nangles), sizeof(*(tr_h)));
204 snew(tr_f, nframes)(tr_f) = save_calloc("tr_f", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 204, (nframes), sizeof(*(tr_f)));
205
206 /* dih[i][j] is the dihedral angle i in frame j */
207 ntrans = 0;
208 for (i = 0; (i < nangles); i++)
209 {
210
211 /*#define OLDIE*/
212#ifdef OLDIE
213 mind = maxd = prev(1-cur) = dih[i][0];
214#else
215 cur_bin = calc_bin(dih[i][0], multiplicity[i], core_frac);
216 rot_occ[cur_bin][i]++;
217#endif
218 for (j = 1; (j < nframes); j++)
219 {
220 new_bin = calc_bin(dih[i][j], multiplicity[i], core_frac);
221 rot_occ[new_bin][i]++;
222#ifndef OLDIE
223 if (cur_bin == 0)
224 {
225 cur_bin = new_bin;
226 }
227 else if ((new_bin != 0) && (cur_bin != new_bin))
228 {
229 cur_bin = new_bin;
230 tr_f[j]++;
231 tr_h[i]++;
232 ntrans++;
233 }
234#else
235 /* why is all this md rubbish periodic? Remove 360 degree periodicity */
236 if ( (dih[i][j] - prev(1-cur)) > M_PI3.14159265358979323846)
237 {
238 dih[i][j] -= 2*M_PI3.14159265358979323846;
239 }
240 else if ( (dih[i][j] - prev(1-cur)) < -M_PI3.14159265358979323846)
241 {
242 dih[i][j] += 2*M_PI3.14159265358979323846;
243 }
244
245 prev(1-cur) = dih[i][j];
246
247 mind = min(mind, dih[i][j])(((mind) < (dih[i][j])) ? (mind) : (dih[i][j]) );
248 maxd = max(maxd, dih[i][j])(((maxd) > (dih[i][j])) ? (maxd) : (dih[i][j]) );
249 if ( (maxd - mind) > 2*M_PI3.14159265358979323846/3) /* or 120 degrees, assuming */
250 { /* multiplicity 3. Not so general.*/
251 tr_f[j]++;
252 tr_h[i]++;
253 maxd = mind = dih[i][j]; /* get ready for next transition */
254 ntrans++;
255 }
256#endif
257 } /* end j */
258 for (k = 0; k < NROT4; k++)
259 {
260 rot_occ[k][i] /= nframes;
261 }
262 } /* end i */
263 fprintf(stderrstderr, "Total number of transitions: %10d\n", ntrans);
264 if (ntrans > 0)
265 {
266 ttime = (dt*nframes*nangles)/ntrans;
267 fprintf(stderrstderr, "Time between transitions: %10.3f ps\n", ttime);
268 }
269
270 /* new by grs - copy transitions from tr_h[] to dlist->ntr[]
271 * and rotamer populations from rot_occ to dlist->rot_occ[]
272 * based on fn histogramming in g_chi. diff roles for i and j here */
273
274 j = 0;
275 for (Dih = 0; (Dih < NONCHI3+maxchi); Dih++)
276 {
277 for (i = 0; (i < nlist); i++)
278 {
279 if (((Dih < edOmega) ) ||
280 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
281 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI3+3] != -1)))
282 {
283 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
284 dlist[i].ntr[Dih] = tr_h[j];
285 for (k = 0; k < NROT4; k++)
286 {
287 dlist[i].rot_occ[Dih][k] = rot_occ[k][j];
288 }
289 j++;
290 }
291 }
292 }
293
294 /* end addition by grs */
295
296 if (bTrans)
297 {
298 sprintf(title, "Number of transitions: %s", grpname);
299 fp = xvgropen(fn_trans, title, "Time (ps)", "# transitions/timeframe", oenv);
300 for (j = 0; (j < nframes); j++)
301 {
302 fprintf(fp, "%10.3f %10d\n", time[j], tr_f[j]);
303 }
304 gmx_ffclose(fp);
305 }
306
307 /* Compute histogram from # transitions per dihedral */
308 /* Use old array */
309 for (j = 0; (j < nframes); j++)
310 {
311 tr_f[j] = 0;
312 }
313 for (i = 0; (i < nangles); i++)
314 {
315 tr_f[tr_h[i]]++;
316 }
317 for (j = nframes; ((tr_f[j-1] == 0) && (j > 0)); j--)
318 {
319 ;
320 }
321
322 ttime = dt*nframes;
323 if (bHisto)
324 {
325 sprintf(title, "Transition time: %s", grpname);
326 fp = xvgropen(fn_histo, title, "Time (ps)", "#", oenv);
327 for (i = j-1; (i > 0); i--)
328 {
329 if (tr_f[i] != 0)
330 {
331 fprintf(fp, "%10.3f %10d\n", ttime/i, tr_f[i]);
332 }
333 }
334 gmx_ffclose(fp);
335 }
336
337 sfree(tr_f)save_free("tr_f", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 337, (tr_f));
338 sfree(tr_h)save_free("tr_h", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 338, (tr_h));
339 for (k = 0; k < NROT4; k++)
340 {
341 sfree(rot_occ[k])save_free("rot_occ[k]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 341, (rot_occ[k]));
342 }
343
344}
345
346void mk_multiplicity_lookup (int *multiplicity, int maxchi,
347 int nlist, t_dlist dlist[], int nangles)
348{
349 /* new by grs - for dihedral j (as in dih[j]) get multiplicity from dlist
350 * and store in multiplicity[j]
351 */
352
353 int j, Dih, i;
354 char name[4];
355
356 j = 0;
357 for (Dih = 0; (Dih < NONCHI3+maxchi); Dih++)
358 {
359 for (i = 0; (i < nlist); i++)
360 {
361 strncpy(name, dlist[i].name, 3)__builtin_strncpy (name, dlist[i].name, 3);
362 name[3] = '\0';
363 if (((Dih < edOmega) ) ||
364 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
365 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI3+3] != -1)))
366 {
367 /* default - we will correct the rest below */
368 multiplicity[j] = 3;
369
370 /* make omegas 2fold, though doesn't make much more sense than 3 */
371 if (Dih == edOmega && (has_dihedral(edOmega, &(dlist[i]))))
372 {
373 multiplicity[j] = 2;
374 }
375
376 /* dihedrals to aromatic rings, COO, CONH2 or guanidinium are 2fold*/
377 if (Dih > edOmega && (dlist[i].atm.Cn[Dih-NONCHI3+3] != -1))
378 {
379 if ( ((strstr(name, "PHE") != NULL((void*)0)) && (Dih == edChi2)) ||
380 ((strstr(name, "TYR") != NULL((void*)0)) && (Dih == edChi2)) ||
381 ((strstr(name, "PTR") != NULL((void*)0)) && (Dih == edChi2)) ||
382 ((strstr(name, "TRP") != NULL((void*)0)) && (Dih == edChi2)) ||
383 ((strstr(name, "HIS") != NULL((void*)0)) && (Dih == edChi2)) ||
384 ((strstr(name, "GLU") != NULL((void*)0)) && (Dih == edChi3)) ||
385 ((strstr(name, "ASP") != NULL((void*)0)) && (Dih == edChi2)) ||
386 ((strstr(name, "GLN") != NULL((void*)0)) && (Dih == edChi3)) ||
387 ((strstr(name, "ASN") != NULL((void*)0)) && (Dih == edChi2)) ||
388 ((strstr(name, "ARG") != NULL((void*)0)) && (Dih == edChi4)) )
389 {
390 multiplicity[j] = 2;
391 }
392 }
393 j++;
394 }
395 }
396 }
397 if (j < nangles)
398 {
399 fprintf(stderrstderr, "WARNING: not all dihedrals found in topology (only %d out of %d)!\n",
400 j, nangles);
401 }
402 /* Check for remaining dihedrals */
403 for (; (j < nangles); j++)
404 {
405 multiplicity[j] = 3;
406 }
407
408}
409
410void mk_chi_lookup (int **lookup, int maxchi,
411 int nlist, t_dlist dlist[])
412{
413
414 /* by grs. should rewrite everything to use this. (but haven't,
415 * and at mmt only used in get_chi_product_traj
416 * returns the dihed number given the residue number (from-0)
417 * and chi (from-0) nr. -1 for chi undefined for that res (eg gly, ala..)*/
418
419 int i, j, Dih, Chi;
420
421 j = 0;
422 for (Dih = 0; (Dih < NONCHI3+maxchi); Dih++)
423 {
424 for (i = 0; (i < nlist); i++)
425 {
426 Chi = Dih - NONCHI3;
427 if (((Dih < edOmega) ) ||
428 ((Dih == edOmega) && (has_dihedral(edOmega, &(dlist[i])))) ||
429 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI3+3] != -1)))
430 {
431 /* grs debug printf("Not OK? i %d j %d Dih %d \n", i, j, Dih) ; */
432 if (Dih > edOmega)
433 {
434 lookup[i][Chi] = j;
435 }
436 j++;
437 }
438 else
439 {
440 lookup[i][Chi] = -1;
441 }
442 }
443 }
444
445}
446
447
448void get_chi_product_traj (real **dih, int nframes, int nlist,
449 int maxchi, t_dlist dlist[], real time[],
450 int **lookup, int *multiplicity, gmx_bool bRb, gmx_bool bNormalize,
451 real core_frac, gmx_bool bAll, const char *fnall,
452 const output_env_t oenv)
453{
454
455 gmx_bool bRotZero, bHaveChi = FALSE0;
456 int accum = 0, index, i, j, k, Xi, n, b;
457 real *chi_prtrj;
458 int *chi_prhist;
459 int nbin;
460 FILE *fp, *fpall;
461 char hisfile[256], histitle[256], *namept;
462
463 int (*calc_bin)(real, int, real);
464
465 /* Analysis of dihedral transitions */
466 fprintf(stderrstderr, "Now calculating Chi product trajectories...\n");
467
468 if (bRb)
469 {
470 calc_bin = calc_RBbin;
471 }
472 else
473 {
474 calc_bin = calc_Nbin;
475 }
476
477 snew(chi_prtrj, nframes)(chi_prtrj) = save_calloc("chi_prtrj", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 477, (nframes), sizeof(*(chi_prtrj)));
478
479 /* file for info on all residues */
480 if (bNormalize)
481 {
482 fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "Probability", oenv);
483 }
484 else
485 {
486 fpall = xvgropen(fnall, "Cumulative Rotamers", "Residue", "# Counts", oenv);
487 }
488
489 for (i = 0; (i < nlist); i++)
490 {
491
492 /* get nbin, the nr. of cumulative rotamers that need to be considered */
493 nbin = 1;
494 for (Xi = 0; Xi < maxchi; Xi++)
495 {
496 index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
497 if (index >= 0)
498 {
499 n = multiplicity[index];
500 nbin = n*nbin;
501 }
502 }
503 nbin += 1; /* for the "zero rotamer", outside the core region */
504
505 for (j = 0; (j < nframes); j++)
506 {
507
508 bRotZero = FALSE0;
509 bHaveChi = TRUE1;
510 index = lookup[i][0]; /* index into dih of chi1 of res i */
511 if (index == -1)
512 {
513 b = 0;
Value stored to 'b' is never read
514 bRotZero = TRUE1;
515 bHaveChi = FALSE0;
516 }
517 else
518 {
519 b = calc_bin(dih[index][j], multiplicity[index], core_frac);
520 accum = b - 1;
521 if (b == 0)
522 {
523 bRotZero = TRUE1;
524 }
525 for (Xi = 1; Xi < maxchi; Xi++)
526 {
527 index = lookup[i][Xi]; /* chi_(Xi+1) of res i (-1 if off end) */
528 if (index >= 0)
529 {
530 n = multiplicity[index];
531 b = calc_bin(dih[index][j], n, core_frac);
532 accum = n * accum + b - 1;
533 if (b == 0)
534 {
535 bRotZero = TRUE1;
536 }
537 }
538 }
539 accum++;
540 }
541 if (bRotZero)
542 {
543 chi_prtrj[j] = 0.0;
544 }
545 else
546 {
547 chi_prtrj[j] = accum;
548 if (accum+1 > nbin)
549 {
550 nbin = accum+1;
551 }
552 }
553 }
554 if (bHaveChi)
555 {
556
557 if (bAll)
558 {
559 /* print cuml rotamer vs time */
560 print_one(oenv, "chiproduct", dlist[i].name, "chi product for",
561 "cumulative rotamer", nframes, time, chi_prtrj);
562 }
563
564 /* make a histogram pf culm. rotamer occupancy too */
565 snew(chi_prhist, nbin)(chi_prhist) = save_calloc("chi_prhist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 565, (nbin), sizeof(*(chi_prhist)));
566 make_histo(NULL((void*)0), nframes, chi_prtrj, nbin, chi_prhist, 0, nbin);
567 if (bAll)
568 {
569 sprintf(hisfile, "histo-chiprod%s.xvg", dlist[i].name);
570 sprintf(histitle, "cumulative rotamer distribution for %s", dlist[i].name);
571 fprintf(stderrstderr, " and %s ", hisfile);
572 fp = xvgropen(hisfile, histitle, "number", "", oenv);
573 fprintf(fp, "@ xaxis tick on\n");
574 fprintf(fp, "@ xaxis tick major 1\n");
575 fprintf(fp, "@ type xy\n");
576 for (k = 0; (k < nbin); k++)
577 {
578 if (bNormalize)
579 {
580 fprintf(fp, "%5d %10g\n", k, (1.0*chi_prhist[k])/nframes);
581 }
582 else
583 {
584 fprintf(fp, "%5d %10d\n", k, chi_prhist[k]);
585 }
586 }
587 fprintf(fp, "&\n");
588 gmx_ffclose(fp);
589 }
590
591 /* and finally print out occupancies to a single file */
592 /* get the gmx from-1 res nr by setting a ptr to the number part
593 * of dlist[i].name - potential bug for 4-letter res names... */
594 namept = dlist[i].name + 3;
595 fprintf(fpall, "%5s ", namept);
596 for (k = 0; (k < nbin); k++)
597 {
598 if (bNormalize)
599 {
600 fprintf(fpall, " %10g", (1.0*chi_prhist[k])/nframes);
601 }
602 else
603 {
604 fprintf(fpall, " %10d", chi_prhist[k]);
605 }
606 }
607 fprintf(fpall, "\n");
608
609 sfree(chi_prhist)save_free("chi_prhist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 609, (chi_prhist));
610 /* histogram done */
611 }
612 }
613
614 sfree(chi_prtrj)save_free("chi_prtrj", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 614, (chi_prtrj));
615 gmx_ffclose(fpall);
616 fprintf(stderrstderr, "\n");
617
618}
619
620void calc_distribution_props(int nh, int histo[], real start,
621 int nkkk, t_karplus kkk[],
622 real *S2)
623{
624 real d, dc, ds, c1, c2, tdc, tds;
625 real fac, ang, invth, Jc;
626 int i, j, th;
627
628 if (nh == 0)
629 {
630 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 630, "No points in histogram (%s, %d)", __FILE__"/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c", __LINE__630);
631 }
632 fac = 2*M_PI3.14159265358979323846/nh;
633
634 /* Compute normalisation factor */
635 th = 0;
636 for (j = 0; (j < nh); j++)
637 {
638 th += histo[j];
639 }
640 invth = 1.0/th;
641
642 for (i = 0; (i < nkkk); i++)
643 {
644 kkk[i].Jc = 0;
645 kkk[i].Jcsig = 0;
646 }
647 tdc = 0, tds = 0;
648 for (j = 0; (j < nh); j++)
649 {
650 d = invth*histo[j];
651 ang = j*fac-start;
652 c1 = cos(ang);
653 c2 = c1*c1;
654 dc = d*c1;
655 ds = d*sin(ang);
656 tdc += dc;
657 tds += ds;
658 for (i = 0; (i < nkkk); i++)
659 {
660 c1 = cos(ang+kkk[i].offset);
661 c2 = c1*c1;
662 Jc = (kkk[i].A*c2 + kkk[i].B*c1 + kkk[i].C);
663 kkk[i].Jc += histo[j]*Jc;
664 kkk[i].Jcsig += histo[j]*sqr(Jc);
665 }
666 }
667 for (i = 0; (i < nkkk); i++)
668 {
669 kkk[i].Jc /= th;
670 kkk[i].Jcsig = sqrt(kkk[i].Jcsig/th-sqr(kkk[i].Jc));
671 }
672 *S2 = tdc*tdc+tds*tds;
673}
674
675static void calc_angles(t_pbc *pbc,
676 int n3, atom_id index[], real ang[], rvec x_s[])
677{
678 int i, ix, t1, t2;
679 rvec r_ij, r_kj;
680 real costh;
681
682 for (i = ix = 0; (ix < n3); i++, ix += 3)
683 {
684 ang[i] = bond_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
685 pbc, r_ij, r_kj, &costh, &t1, &t2);
686 }
687 if (debug)
688 {
689 fprintf(debug, "Angle[0]=%g, costh=%g, index0 = %d, %d, %d\n",
690 ang[0], costh, index[0], index[1], index[2]);
691 pr_rvec(debug, 0, "rij", r_ij, DIM3, TRUE1);
692 pr_rvec(debug, 0, "rkj", r_kj, DIM3, TRUE1);
693 }
694}
695
696static real calc_fraction(real angles[], int nangles)
697{
698 int i;
699 real trans = 0, gauche = 0;
700 real angle;
701
702 for (i = 0; i < nangles; i++)
703 {
704 angle = angles[i] * RAD2DEG(180.0/3.14159265358979323846);
705
706 if (angle > 135 && angle < 225)
707 {
708 trans += 1.0;
709 }
710 else if (angle > 270 && angle < 330)
711 {
712 gauche += 1.0;
713 }
714 else if (angle < 90 && angle > 30)
715 {
716 gauche += 1.0;
717 }
718 }
719 if (trans+gauche > 0)
720 {
721 return trans/(trans+gauche);
722 }
723 else
724 {
725 return 0;
726 }
727}
728
729static void calc_dihs(t_pbc *pbc,
730 int n4, atom_id index[], real ang[], rvec x_s[])
731{
732 int i, ix, t1, t2, t3;
733 rvec r_ij, r_kj, r_kl, m, n;
734 real sign, aaa;
735
736 for (i = ix = 0; (ix < n4); i++, ix += 4)
737 {
738 aaa = dih_angle(x_s[index[ix]], x_s[index[ix+1]], x_s[index[ix+2]],
739 x_s[index[ix+3]], pbc,
740 r_ij, r_kj, r_kl, m, n,
741 &sign, &t1, &t2, &t3);
742
743 ang[i] = aaa; /* not taking into account ryckaert bellemans yet */
744 }
745}
746
747void make_histo(FILE *log,
748 int ndata, real data[], int npoints, int histo[],
749 real minx, real maxx)
750{
751 double dx;
752 int i, ind;
753
754 if (minx == maxx)
755 {
756 minx = maxx = data[0];
757 for (i = 1; (i < ndata); i++)
758 {
759 minx = min(minx, data[i])(((minx) < (data[i])) ? (minx) : (data[i]) );
760 maxx = max(maxx, data[i])(((maxx) > (data[i])) ? (maxx) : (data[i]) );
761 }
762 fprintf(log, "Min data: %10g Max data: %10g\n", minx, maxx);
763 }
764 dx = (double)npoints/(maxx-minx);
765 if (debug)
766 {
767 fprintf(debug,
768 "Histogramming: ndata=%d, nhisto=%d, minx=%g,maxx=%g,dx=%g\n",
769 ndata, npoints, minx, maxx, dx);
770 }
771 for (i = 0; (i < ndata); i++)
772 {
773 ind = (data[i]-minx)*dx;
774 if ((ind >= 0) && (ind < npoints))
775 {
776 histo[ind]++;
777 }
778 else
779 {
780 fprintf(log, "index = %d, data[%d] = %g\n", ind, i, data[i]);
781 }
782 }
783}
784
785void normalize_histo(int npoints, int histo[], real dx, real normhisto[])
786{
787 int i;
788 double d, fac;
789
790 d = 0;
791 for (i = 0; (i < npoints); i++)
792 {
793 d += dx*histo[i];
794 }
795 if (d == 0)
796 {
797 fprintf(stderrstderr, "Empty histogram!\n");
798 return;
799 }
800 fac = 1.0/d;
801 for (i = 0; (i < npoints); i++)
802 {
803 normhisto[i] = fac*histo[i];
804 }
805}
806
807void read_ang_dih(const char *trj_fn,
808 gmx_bool bAngles, gmx_bool bSaveAll, gmx_bool bRb, gmx_bool bPBC,
809 int maxangstat, int angstat[],
810 int *nframes, real **time,
811 int isize, atom_id index[],
812 real **trans_frac,
813 real **aver_angle,
814 real *dih[],
815 const output_env_t oenv)
816{
817 t_pbc *pbc;
818 t_trxstatus *status;
819 int i, angind, natoms, total, teller;
820 int nangles, n_alloc;
821 real t, fraction, pifac, aa, angle;
822 real *angles[2];
823 matrix box;
824 rvec *x;
825 int cur = 0;
826#define prev(1-cur) (1-cur)
827
828 snew(pbc, 1)(pbc) = save_calloc("pbc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 828, (1), sizeof(*(pbc)));
829 natoms = read_first_x(oenv, &status, trj_fn, &t, &x, box);
830
831 if (bAngles)
832 {
833 nangles = isize/3;
834 pifac = M_PI3.14159265358979323846;
835 }
836 else
837 {
838 nangles = isize/4;
839 pifac = 2.0*M_PI3.14159265358979323846;
840 }
841 snew(angles[cur], nangles)(angles[cur]) = save_calloc("angles[cur]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 841, (nangles), sizeof(*(angles[cur])));
842 snew(angles[prev], nangles)(angles[(1-cur)]) = save_calloc("angles[prev]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 842, (nangles), sizeof(*(angles[(1-cur)])));
843
844 /* Start the loop over frames */
845 total = 0;
846 teller = 0;
847 n_alloc = 0;
848 *time = NULL((void*)0);
849 *trans_frac = NULL((void*)0);
850 *aver_angle = NULL((void*)0);
851
852 do
853 {
854 if (teller >= n_alloc)
855 {
856 n_alloc += 100;
857 if (bSaveAll)
858 {
859 for (i = 0; (i < nangles); i++)
860 {
861 srenew(dih[i], n_alloc)(dih[i]) = save_realloc("dih[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 861, (dih[i]), (n_alloc), sizeof(*(dih[i])));
862 }
863 }
864 srenew(*time, n_alloc)(*time) = save_realloc("*time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 864, (*time), (n_alloc), sizeof(*(*time)));
865 srenew(*trans_frac, n_alloc)(*trans_frac) = save_realloc("*trans_frac", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 865, (*trans_frac), (n_alloc), sizeof(*(*trans_frac)));
866 srenew(*aver_angle, n_alloc)(*aver_angle) = save_realloc("*aver_angle", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 866, (*aver_angle), (n_alloc), sizeof(*(*aver_angle)));
867 }
868
869 (*time)[teller] = t;
870
871 if (pbc)
872 {
873 set_pbc(pbc, -1, box);
874 }
875
876 if (bAngles)
877 {
878 calc_angles(pbc, isize, index, angles[cur], x);
879 }
880 else
881 {
882 calc_dihs(pbc, isize, index, angles[cur], x);
883
884 /* Trans fraction */
885 fraction = calc_fraction(angles[cur], nangles);
886 (*trans_frac)[teller] = fraction;
887
888 /* Change Ryckaert-Bellemans dihedrals to polymer convention
889 * Modified 990913 by Erik:
890 * We actually shouldn't change the convention, since it's
891 * calculated from polymer above, but we change the intervall
892 * from [-180,180] to [0,360].
893 */
894 if (bRb)
895 {
896 for (i = 0; (i < nangles); i++)
897 {
898 if (angles[cur][i] <= 0.0)
899 {
900 angles[cur][i] += 2*M_PI3.14159265358979323846;
901 }
902 }
903 }
904
905 /* Periodicity in dihedral space... */
906 if (bPBC)
907 {
908 for (i = 0; (i < nangles); i++)
909 {
910 real dd = angles[cur][i];
911 angles[cur][i] = atan2(sin(dd), cos(dd));
912 }
913 }
914 else
915 {
916 if (teller > 1)
917 {
918 for (i = 0; (i < nangles); i++)
919 {
920 while (angles[cur][i] <= angles[prev(1-cur)][i] - M_PI3.14159265358979323846)
921 {
922 angles[cur][i] += 2*M_PI3.14159265358979323846;
923 }
924 while (angles[cur][i] > angles[prev(1-cur)][i] + M_PI3.14159265358979323846)
925 {
926 angles[cur][i] -= 2*M_PI3.14159265358979323846;
927 }
928 }
929 }
930 }
931 }
932
933 /* Average angles */
934 aa = 0;
935 for (i = 0; (i < nangles); i++)
936 {
937 aa = aa+angles[cur][i];
938
939 /* angle in rad / 2Pi * max determines bin. bins go from 0 to maxangstat,
940 even though scale goes from -pi to pi (dihedral) or -pi/2 to pi/2
941 (angle) Basically: translate the x-axis by Pi. Translate it back by
942 -Pi when plotting.
943 */
944
945 angle = angles[cur][i];
946 if (!bAngles)
947 {
948 while (angle < -M_PI3.14159265358979323846)
949 {
950 angle += 2*M_PI3.14159265358979323846;
951 }
952 while (angle >= M_PI3.14159265358979323846)
953 {
954 angle -= 2*M_PI3.14159265358979323846;
955 }
956
957 angle += M_PI3.14159265358979323846;
958 }
959
960 /* Update the distribution histogram */
961 angind = (int) ((angle*maxangstat)/pifac + 0.5);
962 if (angind == maxangstat)
963 {
964 angind = 0;
965 }
966 if ( (angind < 0) || (angind >= maxangstat) )
967 {
968 /* this will never happen */
969 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 969, "angle (%f) index out of range (0..%d) : %d\n",
970 angle, maxangstat, angind);
971 }
972
973 angstat[angind]++;
974 if (angind == maxangstat)
975 {
976 fprintf(stderrstderr, "angle %d fr %d = %g\n", i, cur, angle);
977 }
978
979 total++;
980 }
981
982 /* average over all angles */
983 (*aver_angle)[teller] = (aa/nangles);
984
985 /* this copies all current dih. angles to dih[i], teller is frame */
986 if (bSaveAll)
987 {
988 for (i = 0; i < nangles; i++)
989 {
990 dih[i][teller] = angles[cur][i];
991 }
992 }
993
994 /* Swap buffers */
995 cur = prev(1-cur);
996
997 /* Increment loop counter */
998 teller++;
999 }
1000 while (read_next_x(oenv, status, &t, x, box));
1001 close_trj(status);
1002
1003 sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 1003, (x));
1004 sfree(angles[cur])save_free("angles[cur]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 1004, (angles[cur]));
1005 sfree(angles[prev])save_free("angles[prev]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/anadih.c"
, 1005, (angles[(1-cur)]));
1006
1007 *nframes = teller;
1008}