File: | gromacs/gmxana/anadih.c |
Location: | line 689, column 9 |
Description: | Function call argument is an uninitialized value |
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; | |||
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 | } |