File: | gromacs/gmxana/anadih.c |
Location: | line 653, column 9 |
Description: | Value stored to 'c2' 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 | #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 | |
58 | void 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 | |
77 | static 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 | |
100 | static 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 | |
132 | void 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 | |
159 | void 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 | |
346 | void 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 | |
410 | void 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 | |
448 | void 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; |
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 | |
620 | void 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; |
Value stored to 'c2' is never read | |
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 | |
675 | static 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 | |
696 | static 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 | |
729 | static 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 | |
747 | void 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 | |
785 | void 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 | |
807 | void 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 | } |