Bug Summary

File:gromacs/gmxana/gmx_angle.c
Location:line 194, column 13
Description:Array access results in a null pointer dereference

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#include <math.h>
41#include <string.h>
42
43#include "physics.h"
44#include "typedefs.h"
45#include "gromacs/utility/smalloc.h"
46#include "gromacs/utility/futil.h"
47#include "gromacs/commandline/pargs.h"
48#include "copyrite.h"
49#include "gromacs/math/vec.h"
50#include "index.h"
51#include "macros.h"
52#include "gromacs/utility/fatalerror.h"
53#include "gromacs/fileio/xvgr.h"
54#include "viewit.h"
55#include "gstat.h"
56#include "gromacs/fileio/trnio.h"
57#include "gmx_ana.h"
58
59
60static void dump_dih_trn(int nframes, int nangles, real **dih, const char *fn,
61 real *time)
62{
63 int i, j, k, l, m, na;
64 t_fileio *trn;
65 rvec *x;
66 matrix box = {{2, 0, 0}, {0, 2, 0}, {0, 0, 2}};
67
68 na = (nangles*2);
69 if ((na % 3) != 0)
70 {
71 na = 1+na/3;
72 }
73 else
74 {
75 na = na/3;
76 }
77 printf("There are %d dihedrals. Will fill %d atom positions with cos/sin\n",
78 nangles, na);
79 snew(x, na)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 79, (na), sizeof(*(x)));
80 trn = open_trn(fn, "w");
81 for (i = 0; (i < nframes); i++)
82 {
83 k = l = 0;
84 for (j = 0; (j < nangles); j++)
85 {
86 for (m = 0; (m < 2); m++)
87 {
88 x[k][l] = (m == 0) ? cos(dih[j][i]) : sin(dih[j][i]);
89 l++;
90 if (l == DIM3)
91 {
92 l = 0;
93 k++;
94 }
95 }
96 }
97 fwrite_trn(trn, i, time[i], 0, box, na, x, NULL((void*)0), NULL((void*)0));
98 }
99 close_trn(trn);
100 sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 100, (x));
101}
102
103int gmx_g_angle(int argc, char *argv[])
104{
105 static const char *desc[] = {
106 "[THISMODULE] computes the angle distribution for a number of angles",
107 "or dihedrals.[PAR]",
108 "With option [TT]-ov[tt], you can plot the average angle of",
109 "a group of angles as a function of time. With the [TT]-all[tt] option,",
110 "the first graph is the average and the rest are the individual angles.[PAR]",
111 "With the [TT]-of[tt] option, [THISMODULE] also calculates the fraction of trans",
112 "dihedrals (only for dihedrals) as function of time, but this is",
113 "probably only fun for a select few.[PAR]",
114 "With option [TT]-oc[tt], a dihedral correlation function is calculated.[PAR]",
115 "It should be noted that the index file must contain",
116 "atom triplets for angles or atom quadruplets for dihedrals.",
117 "If this is not the case, the program will crash.[PAR]",
118 "With option [TT]-or[tt], a trajectory file is dumped containing cos and",
119 "sin of selected dihedral angles, which subsequently can be used as",
120 "input for a principal components analysis using [gmx-covar].[PAR]",
121 "Option [TT]-ot[tt] plots when transitions occur between",
122 "dihedral rotamers of multiplicity 3 and [TT]-oh[tt]",
123 "records a histogram of the times between such transitions,",
124 "assuming the input trajectory frames are equally spaced in time."
125 };
126 static const char *opt[] = { NULL((void*)0), "angle", "dihedral", "improper", "ryckaert-bellemans", NULL((void*)0) };
127 static gmx_bool bALL = FALSE0, bChandler = FALSE0, bAverCorr = FALSE0, bPBC = TRUE1;
128 static real binwidth = 1;
129 t_pargs pa[] = {
130 { "-type", FALSE0, etENUM, {opt},
131 "Type of angle to analyse" },
132 { "-all", FALSE0, etBOOL, {&bALL},
133 "Plot all angles separately in the averages file, in the order of appearance in the index file." },
134 { "-binwidth", FALSE0, etREAL, {&binwidth},
135 "binwidth (degrees) for calculating the distribution" },
136 { "-periodic", FALSE0, etBOOL, {&bPBC},
137 "Print dihedral angles modulo 360 degrees" },
138 { "-chandler", FALSE0, etBOOL, {&bChandler},
139 "Use Chandler correlation function (N[trans] = 1, N[gauche] = 0) rather than cosine correlation function. Trans is defined as phi < -60 or phi > 60." },
140 { "-avercorr", FALSE0, etBOOL, {&bAverCorr},
141 "Average the correlation functions for the individual angles/dihedrals" }
142 };
143 static const char *bugs[] = {
144 "Counting transitions only works for dihedrals with multiplicity 3"
145 };
146
147 FILE *out;
148 real tmp, dt;
149 int status, isize;
150 atom_id *index;
151 char *grpname;
152 real maxang, Jc, S2, norm_fac, maxstat;
153 unsigned long mode;
154 int nframes, maxangstat, mult, *angstat;
155 int i, j, total, nangles, natoms, nat2, first, last, angind;
156 gmx_bool bAver, bRb, bPeriodic,
157 bFrac, /* calculate fraction too? */
158 bTrans, /* worry about transtions too? */
159 bCorr; /* correlation function ? */
160 real t, aa, aver, aver2, aversig, fraction; /* fraction trans dihedrals */
161 double tfrac = 0;
162 char title[256];
163 real **dih = NULL((void*)0); /* mega array with all dih. angles at all times*/
164 char buf[80];
165 real *time, *trans_frac, *aver_angle;
166 t_filenm fnm[] = {
167 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
168 { efNDX, NULL((void*)0), "angle", ffREAD1<<1 },
169 { efXVG, "-od", "angdist", ffWRITE1<<2 },
170 { efXVG, "-ov", "angaver", ffOPTWR(1<<2| 1<<3) },
171 { efXVG, "-of", "dihfrac", ffOPTWR(1<<2| 1<<3) },
172 { efXVG, "-ot", "dihtrans", ffOPTWR(1<<2| 1<<3) },
173 { efXVG, "-oh", "trhisto", ffOPTWR(1<<2| 1<<3) },
174 { efXVG, "-oc", "dihcorr", ffOPTWR(1<<2| 1<<3) },
175 { efTRR, "-or", NULL((void*)0), ffOPTWR(1<<2| 1<<3) }
176 };
177#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
178 int npargs;
179 t_pargs *ppa;
180 output_env_t oenv;
181
182 npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])));
183 ppa = add_acf_pargs(&npargs, pa);
184 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW(1<<5) | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_BE_NICE(1<<13),
1
Taking false branch
→
185 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, ppa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, asize(bugs)((int)(sizeof(bugs)/sizeof((bugs)[0]))), bugs,
186 &oenv))
187 {
188 return 0;
189 }
190
191 mult = 4;
192 maxang = 360.0;
193 bRb = FALSE0;
194 switch (opt[0][0])
2
←
Array access results in a null pointer dereference
195 {
196 case 'a':
197 mult = 3;
198 maxang = 180.0;
199 break;
200 case 'd':
201 break;
202 case 'i':
203 break;
204 case 'r':
205 bRb = TRUE1;
206 break;
207 }
208
209 if (opt2bSet("-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
210 {
211 if (mult != 4)
212 {
213 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 213, "Can not combine angles with trn dump");
214 }
215 else
216 {
217 please_cite(stdoutstdout, "Mu2005a");
218 }
219 }
220
221 /* Calculate bin size */
222 maxangstat = (int)(maxang/binwidth+0.5);
223 binwidth = maxang/maxangstat;
224
225 rd_index(ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &isize, &index, &grpname);
226 nangles = isize/mult;
227 if ((isize % mult) != 0)
228 {
229 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 229, "number of index elements not multiple of %d, "
230 "these can not be %s\n",
231 mult, (mult == 3) ? "angle triplets" : "dihedral quadruplets");
232 }
233
234
235 /* Check whether specific analysis has to be performed */
236 bCorr = opt2bSet("-oc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
237 bAver = opt2bSet("-ov", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
238 bTrans = opt2bSet("-ot", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
239 bFrac = opt2bSet("-of", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
240 if (bTrans && opt[0][0] != 'd')
241 {
242 fprintf(stderrstderr, "Option -ot should only accompany -type dihedral. Disabling -ot.\n");
243 bTrans = FALSE0;
244 }
245
246 if (bChandler && !bCorr)
247 {
248 bCorr = TRUE1;
249 }
250
251 if (bFrac && !bRb)
252 {
253 fprintf(stderrstderr, "Warning:"
254 " calculating fractions as defined in this program\n"
255 "makes sense for Ryckaert Bellemans dihs. only. Ignoring -of\n\n");
256 bFrac = FALSE0;
257 }
258
259 if ( (bTrans || bFrac || bCorr) && mult == 3)
260 {
261 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 261, "Can only do transition, fraction or correlation\n"
262 "on dihedrals. Select -d\n");
263 }
264
265 /*
266 * We need to know the nr of frames so we can allocate memory for an array
267 * with all dihedral angles at all timesteps. Works for me.
268 */
269 if (bTrans || bCorr || bALL || opt2bSet("-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
270 {
271 snew(dih, nangles)(dih) = save_calloc("dih", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 271, (nangles), sizeof(*(dih)));
272 }
273
274 snew(angstat, maxangstat)(angstat) = save_calloc("angstat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 274, (maxangstat), sizeof(*(angstat)));
275
276 read_ang_dih(ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), (mult == 3),
277 bALL || bCorr || bTrans || opt2bSet("-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
278 bRb, bPBC, maxangstat, angstat,
279 &nframes, &time, isize, index, &trans_frac, &aver_angle, dih,
280 oenv);
281
282 dt = (time[nframes-1]-time[0])/(nframes-1);
283
284 if (bAver)
285 {
286 sprintf(title, "Average Angle: %s", grpname);
287 out = xvgropen(opt2fn("-ov", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
288 title, "Time (ps)", "Angle (degrees)", oenv);
289 for (i = 0; (i < nframes); i++)
290 {
291 fprintf(out, "%10.5f %8.3f", time[i], aver_angle[i]*RAD2DEG(180.0/3.14159265358979323846));
292 if (bALL)
293 {
294 for (j = 0; (j < nangles); j++)
295 {
296 if (bPBC)
297 {
298 real dd = dih[j][i];
299 fprintf(out, " %8.3f", atan2(sin(dd), cos(dd))*RAD2DEG(180.0/3.14159265358979323846));
300 }
301 else
302 {
303 fprintf(out, " %8.3f", dih[j][i]*RAD2DEG(180.0/3.14159265358979323846));
304 }
305 }
306 }
307 fprintf(out, "\n");
308 }
309 gmx_ffclose(out);
310 }
311 if (opt2bSet("-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
312 {
313 dump_dih_trn(nframes, nangles, dih, opt2fn("-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), time);
314 }
315
316 if (bFrac)
317 {
318 sprintf(title, "Trans fraction: %s", grpname);
319 out = xvgropen(opt2fn("-of", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
320 title, "Time (ps)", "Fraction", oenv);
321 tfrac = 0.0;
322 for (i = 0; (i < nframes); i++)
323 {
324 fprintf(out, "%10.5f %10.3f\n", time[i], trans_frac[i]);
325 tfrac += trans_frac[i];
326 }
327 gmx_ffclose(out);
328
329 tfrac /= nframes;
330 fprintf(stderrstderr, "Average trans fraction: %g\n", tfrac);
331 }
332 sfree(trans_frac)save_free("trans_frac", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_angle.c"
, 332, (trans_frac));
333
334 if (bTrans)
335 {
336 ana_dih_trans(opt2fn("-ot", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-oh", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
337 dih, nframes, nangles, grpname, time, bRb, oenv);
338 }
339
340 if (bCorr)
341 {
342 /* Autocorrelation function */
343 if (nframes < 2)
344 {
345 fprintf(stderrstderr, "Not enough frames for correlation function\n");
346 }
347 else
348 {
349
350 if (bChandler)
351 {
352 real dval, sixty = DEG2RAD(3.14159265358979323846/180.0)*60;
353 gmx_bool bTest;
354
355 for (i = 0; (i < nangles); i++)
356 {
357 for (j = 0; (j < nframes); j++)
358 {
359 dval = dih[i][j];
360 if (bRb)
361 {
362 bTest = (dval > -sixty) && (dval < sixty);
363 }
364 else
365 {
366 bTest = (dval < -sixty) || (dval > sixty);
367 }
368 if (bTest)
369 {
370 dih[i][j] = dval-tfrac;
371 }
372 else
373 {
374 dih[i][j] = -tfrac;
375 }
376 }
377 }
378 }
379 if (bChandler)
380 {
381 mode = eacNormal(1<<0);
382 }
383 else
384 {
385 mode = eacCos(1<<1);
386 }
387 do_autocorr(opt2fn("-oc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), oenv,
388 "Dihedral Autocorrelation Function",
389 nframes, nangles, dih, dt, mode, bAverCorr);
390 }
391 }
392
393
394 /* Determine the non-zero part of the distribution */
395 for (first = 0; (first < maxangstat-1) && (angstat[first+1] == 0); first++)
396 {
397 ;
398 }
399 for (last = maxangstat-1; (last > 0) && (angstat[last-1] == 0); last--)
400 {
401 ;
402 }
403
404 aver = aver2 = 0;
405 for (i = 0; (i < nframes); i++)
406 {
407 aver += RAD2DEG(180.0/3.14159265358979323846)*aver_angle[i];
408 aver2 += sqr(RAD2DEG(180.0/3.14159265358979323846)*aver_angle[i]);
409 }
410 aver /= (real) nframes;
411 aver2 /= (real) nframes;
412 aversig = sqrt(aver2-sqr(aver));
413 printf("Found points in the range from %d to %d (max %d)\n",
414 first, last, maxangstat);
415 printf(" < angle > = %g\n", aver);
416 printf("< angle^2 > = %g\n", aver2);
417 printf("Std. Dev. = %g\n", aversig);
418
419 if (mult == 3)
420 {
421 sprintf(title, "Angle Distribution: %s", grpname);
422 }
423 else
424 {
425 sprintf(title, "Dihedral Distribution: %s", grpname);
426
427 calc_distribution_props(maxangstat, angstat, -180.0, 0, NULL((void*)0), &S2);
428 fprintf(stderrstderr, "Order parameter S^2 = %g\n", S2);
429 }
430
431 bPeriodic = (mult == 4) && (first == 0) && (last == maxangstat-1);
432
433 out = xvgropen(opt2fn("-od", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, "Degrees", "", oenv);
434 if (output_env_get_print_xvgr_codes(oenv))
435 {
436 fprintf(out, "@ subtitle \"average angle: %g\\So\\N\"\n", aver);
437 }
438 norm_fac = 1.0/(nangles*nframes*binwidth);
439 if (bPeriodic)
440 {
441 maxstat = 0;
442 for (i = first; (i <= last); i++)
443 {
444 maxstat = max(maxstat, angstat[i]*norm_fac)(((maxstat) > (angstat[i]*norm_fac)) ? (maxstat) : (angstat
[i]*norm_fac) );
445 }
446 fprintf(out, "@with g0\n");
447 fprintf(out, "@ world xmin -180\n");
448 fprintf(out, "@ world xmax 180\n");
449 fprintf(out, "@ world ymin 0\n");
450 fprintf(out, "@ world ymax %g\n", maxstat*1.05);
451 fprintf(out, "@ xaxis tick major 60\n");
452 fprintf(out, "@ xaxis tick minor 30\n");
453 fprintf(out, "@ yaxis tick major 0.005\n");
454 fprintf(out, "@ yaxis tick minor 0.0025\n");
455 }
456 for (i = first; (i <= last); i++)
457 {
458 fprintf(out, "%10g %10f\n", i*binwidth+180.0-maxang, angstat[i]*norm_fac);
459 }
460 if (bPeriodic)
461 {
462 /* print first bin again as last one */
463 fprintf(out, "%10g %10f\n", 180.0, angstat[0]*norm_fac);
464 }
465
466 gmx_ffclose(out);
467
468 do_view(oenv, opt2fn("-od", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
469 if (bAver)
470 {
471 do_view(oenv, opt2fn("-ov", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
472 }
473
474 return 0;
475}