Bug Summary

File:gromacs/gmxana/gmx_rdf.c
Location:line 915, column 16
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
38#ifdef HAVE_CONFIG_H1
39#include <config.h>
40#endif
41
42#include <math.h>
43
44#include "typedefs.h"
45#include "macros.h"
46#include "gromacs/math/vec.h"
47#include "pbc.h"
48#include "gromacs/fileio/xvgr.h"
49#include "viewit.h"
50#include "gromacs/utility/futil.h"
51#include "gromacs/commandline/pargs.h"
52#include "gromacs/fileio/tpxio.h"
53#include "gromacs/fileio/trxio.h"
54#include "physics.h"
55#include "index.h"
56#include "gromacs/utility/smalloc.h"
57#include "calcgrid.h"
58#include "nrnb.h"
59#include "coulomb.h"
60#include "gstat.h"
61#include "gromacs/fileio/matio.h"
62#include "gmx_ana.h"
63#include "names.h"
64
65#include "gromacs/utility/fatalerror.h"
66
67static void check_box_c(matrix box)
68{
69 if (fabs(box[ZZ2][XX0]) > GMX_REAL_EPS5.96046448E-08*box[ZZ2][ZZ2] ||
70 fabs(box[ZZ2][YY1]) > GMX_REAL_EPS5.96046448E-08*box[ZZ2][ZZ2])
71 {
72 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 72,
73 "The last box vector is not parallel to the z-axis: %f %f %f",
74 box[ZZ2][XX0], box[ZZ2][YY1], box[ZZ2][ZZ2]);
75 }
76}
77
78static void calc_comg(int is, int *coi, int *index, gmx_bool bMass, t_atom *atom,
79 rvec *x, rvec *x_comg)
80{
81 int c, i, d;
82 rvec xc;
83 real mtot, m;
84
85 if (bMass && atom == NULL((void*)0))
86 {
87 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 87, "No masses available while mass weighting was requested");
88 }
89
90 for (c = 0; c < is; c++)
91 {
92 clear_rvec(xc);
93 mtot = 0;
94 for (i = coi[c]; i < coi[c+1]; i++)
95 {
96 if (bMass)
97 {
98 m = atom[index[i]].m;
99 for (d = 0; d < DIM3; d++)
100 {
101 xc[d] += m*x[index[i]][d];
102 }
103 mtot += m;
104 }
105 else
106 {
107 rvec_inc(xc, x[index[i]]);
108 mtot += 1.0;
109 }
110 }
111 svmul(1/mtot, xc, x_comg[c]);
112 }
113}
114
115static void split_group(int isize, int *index, char *grpname,
116 t_topology *top, char type,
117 int *is_out, int **coi_out)
118{
119 t_block *mols = NULL((void*)0);
120 t_atom *atom = NULL((void*)0);
121 int is, *coi;
122 int cur, mol, res, i, a, i1;
123
124 /* Split up the group in molecules or residues */
125 switch (type)
126 {
127 case 'm':
128 mols = &top->mols;
129 break;
130 case 'r':
131 atom = top->atoms.atom;
132 break;
133 default:
134 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 134, "Unknown rdf option '%s'", type);
135 }
136 snew(coi, isize+1)(coi) = save_calloc("coi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 136, (isize+1), sizeof(*(coi)));
137 is = 0;
138 cur = -1;
139 mol = 0;
140 for (i = 0; i < isize; i++)
141 {
142 a = index[i];
143 if (type == 'm')
144 {
145 /* Check if the molecule number has changed */
146 i1 = mols->index[mol+1];
147 while (a >= i1)
148 {
149 mol++;
150 i1 = mols->index[mol+1];
151 }
152 if (mol != cur)
153 {
154 coi[is++] = i;
155 cur = mol;
156 }
157 }
158 else if (type == 'r')
159 {
160 /* Check if the residue index has changed */
161 res = atom[a].resind;
162 if (res != cur)
163 {
164 coi[is++] = i;
165 cur = res;
166 }
167 }
168 }
169 coi[is] = i;
170 srenew(coi, is+1)(coi) = save_realloc("coi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 170, (coi), (is+1), sizeof(*(coi)));
171 printf("Group '%s' of %d atoms consists of %d %s\n",
172 grpname, isize, is,
173 (type == 'm' ? "molecules" : "residues"));
174
175 *is_out = is;
176 *coi_out = coi;
177}
178
179static void do_rdf(const char *fnNDX, const char *fnTPS, const char *fnTRX,
180 const char *fnRDF, const char *fnCNRDF, const char *fnHQ,
181 gmx_bool bCM, const char *close,
182 const char **rdft, gmx_bool bXY, gmx_bool bPBC, gmx_bool bNormalize,
183 real cutoff, real binwidth, real fade, int ng,
184 const output_env_t oenv)
185{
186 FILE *fp;
187 t_trxstatus *status;
188 char outf1[STRLEN4096], outf2[STRLEN4096];
189 char title[STRLEN4096], gtitle[STRLEN4096], refgt[30];
190 int g, natoms, i, ii, j, k, nbin, j0, j1, n, nframes;
191 int **count;
192 char **grpname;
193 int *isize, isize_cm = 0, nrdf = 0, max_i, isize0, isize_g;
194 atom_id **index, *index_cm = NULL((void*)0);
195 gmx_int64_t *sum;
196 real t, rmax2, cut2, r, r2, r2ii, invhbinw, normfac;
197 real segvol, spherevol, prev_spherevol, **rdf;
198 rvec *x, dx, *x0 = NULL((void*)0), *x_i1, xi;
199 real *inv_segvol, invvol, invvol_sum, rho;
200 gmx_bool bClose, *bExcl, bTop, bNonSelfExcl;
201 matrix box, box_pbc;
202 int **npairs;
203 atom_id ix, jx, ***pairs;
204 t_topology *top = NULL((void*)0);
205 int ePBC = -1, ePBCrdf = -1;
206 t_block *mols = NULL((void*)0);
207 t_blocka *excl;
208 t_atom *atom = NULL((void*)0);
209 t_pbc pbc;
210 gmx_rmpbc_t gpbc = NULL((void*)0);
211 int *is = NULL((void*)0), **coi = NULL((void*)0), cur, mol, i1, res, a;
212
213 excl = NULL((void*)0);
214
215 bClose = (close[0] != 'n');
216
217 if (fnTPS)
218 {
219 snew(top, 1)(top) = save_calloc("top", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 219, (1), sizeof(*(top)));
220 bTop = read_tps_conf(fnTPS, title, top, &ePBC, &x, NULL((void*)0), box, TRUE1);
221 if (bTop && !(bCM || bClose))
222 {
223 /* get exclusions from topology */
224 excl = &(top->excls);
225 }
226 }
227 snew(grpname, ng+1)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 227, (ng+1), sizeof(*(grpname)));
228 snew(isize, ng+1)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 228, (ng+1), sizeof(*(isize)));
229 snew(index, ng+1)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 229, (ng+1), sizeof(*(index)));
230 fprintf(stderrstderr, "\nSelect a reference group and %d group%s\n",
231 ng, ng == 1 ? "" : "s");
232 if (fnTPS)
233 {
234 get_index(&(top->atoms), fnNDX, ng+1, isize, index, grpname);
235 atom = top->atoms.atom;
236 }
237 else
238 {
239 rd_index(fnNDX, ng+1, isize, index, grpname);
240 }
241
242 if (bCM || bClose)
243 {
244 snew(is, 1)(is) = save_calloc("is", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 244, (1), sizeof(*(is)));
245 snew(coi, 1)(coi) = save_calloc("coi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 245, (1), sizeof(*(coi)));
246 if (bClose)
247 {
248 split_group(isize[0], index[0], grpname[0], top, close[0], &is[0], &coi[0]);
249 }
250 }
251 if (rdft[0][0] != 'a')
252 {
253 /* Split up all the groups in molecules or residues */
254 srenew(is, ng+1)(is) = save_realloc("is", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 254, (is), (ng+1), sizeof(*(is)));
255 srenew(coi, ng+1)(coi) = save_realloc("coi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 255, (coi), (ng+1), sizeof(*(coi)));
256 for (g = ((bCM || bClose) ? 1 : 0); g < ng+1; g++)
257 {
258 split_group(isize[g], index[g], grpname[g], top, rdft[0][0], &is[g], &coi[g]);
259 }
260 }
261
262 if (bCM)
263 {
264 is[0] = 1;
265 snew(coi[0], is[0]+1)(coi[0]) = save_calloc("coi[0]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 265, (is[0]+1), sizeof(*(coi[0])));
266 coi[0][0] = 0;
267 coi[0][1] = isize[0];
268 isize0 = is[0];
269 snew(x0, isize0)(x0) = save_calloc("x0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 269, (isize0), sizeof(*(x0)));
270 }
271 else if (bClose || rdft[0][0] != 'a')
272 {
273 isize0 = is[0];
274 snew(x0, isize0)(x0) = save_calloc("x0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 274, (isize0), sizeof(*(x0)));
275 }
276 else
277 {
278 isize0 = isize[0];
279 }
280
281 natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
282 if (!natoms)
283 {
284 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 284, "Could not read coordinates from statusfile\n");
285 }
286 if (fnTPS)
287 {
288 /* check with topology */
289 if (natoms > top->atoms.nr)
290 {
291 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 291, "Trajectory (%d atoms) does not match topology (%d atoms)",
292 natoms, top->atoms.nr);
293 }
294 }
295 /* check with index groups */
296 for (i = 0; i < ng+1; i++)
297 {
298 for (j = 0; j < isize[i]; j++)
299 {
300 if (index[i][j] >= natoms)
301 {
302 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 302, "Atom index (%d) in index group %s (%d atoms) larger "
303 "than number of atoms in trajectory (%d atoms)",
304 index[i][j], grpname[i], isize[i], natoms);
305 }
306 }
307 }
308
309 /* initialize some handy things */
310 if (ePBC == -1)
311 {
312 ePBC = guess_ePBC(box);
313 }
314 copy_mat(box, box_pbc);
315 if (bXY)
316 {
317 check_box_c(box);
318 switch (ePBC)
319 {
320 case epbcXYZ:
321 case epbcXY: ePBCrdf = epbcXY; break;
322 case epbcNONE: ePBCrdf = epbcNONE; break;
323 default:
324 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 324, "xy-rdf's are not supported for pbc type'%s'",
325 EPBC(ePBC)((((ePBC) < 0) || ((ePBC) >= (epbcNR))) ? "UNDEFINED" :
(epbc_names)[ePBC]));
326 break;
327 }
328 /* Make sure the z-height does not influence the cut-off */
329 box_pbc[ZZ2][ZZ2] = 2*max(box[XX][XX], box[YY][YY])(((box[0][0]) > (box[1][1])) ? (box[0][0]) : (box[1][1]) );
330 }
331 else
332 {
333 ePBCrdf = ePBC;
334 }
335 if (bPBC)
336 {
337 rmax2 = 0.99*0.99*max_cutoff2(bXY ? epbcXY : epbcXYZ, box_pbc);
338 }
339 else
340 {
341 rmax2 = sqr(3*max(box[XX][XX], max(box[YY][YY], box[ZZ][ZZ]))(((box[0][0]) > ((((box[1][1]) > (box[2][2])) ? (box[1]
[1]) : (box[2][2]) ))) ? (box[0][0]) : ((((box[1][1]) > (box
[2][2])) ? (box[1][1]) : (box[2][2]) )) ));
342 }
343 if (debug)
344 {
345 fprintf(debug, "rmax2 = %g\n", rmax2);
346 }
347
348 /* We use the double amount of bins, so we can correctly
349 * write the rdf and rdf_cn output at i*binwidth values.
350 */
351 nbin = (int)(sqrt(rmax2) * 2 / binwidth);
352 invhbinw = 2.0 / binwidth;
353 cut2 = sqr(cutoff);
354
355 snew(count, ng)(count) = save_calloc("count", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 355, (ng), sizeof(*(count)));
356 snew(pairs, ng)(pairs) = save_calloc("pairs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 356, (ng), sizeof(*(pairs)));
357 snew(npairs, ng)(npairs) = save_calloc("npairs", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 357, (ng), sizeof(*(npairs)));
358
359 snew(bExcl, natoms)(bExcl) = save_calloc("bExcl", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 359, (natoms), sizeof(*(bExcl)));
360 max_i = 0;
361 for (g = 0; g < ng; g++)
362 {
363 if (isize[g+1] > max_i)
364 {
365 max_i = isize[g+1];
366 }
367
368 /* this is THE array */
369 snew(count[g], nbin+1)(count[g]) = save_calloc("count[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 369, (nbin+1), sizeof(*(count[g])));
370
371 /* make pairlist array for groups and exclusions */
372 snew(pairs[g], isize[0])(pairs[g]) = save_calloc("pairs[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 372, (isize[0]), sizeof(*(pairs[g])));
373 snew(npairs[g], isize[0])(npairs[g]) = save_calloc("npairs[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 373, (isize[0]), sizeof(*(npairs[g])));
374 for (i = 0; i < isize[0]; i++)
375 {
376 /* We can only have exclusions with atomic rdfs */
377 if (!(bCM || bClose || rdft[0][0] != 'a'))
378 {
379 ix = index[0][i];
380 for (j = 0; j < natoms; j++)
381 {
382 bExcl[j] = FALSE0;
383 }
384 /* exclusions? */
385 if (excl)
386 {
387 for (j = excl->index[ix]; j < excl->index[ix+1]; j++)
388 {
389 bExcl[excl->a[j]] = TRUE1;
390 }
391 }
392 k = 0;
393 snew(pairs[g][i], isize[g+1])(pairs[g][i]) = save_calloc("pairs[g][i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 393, (isize[g+1]), sizeof(*(pairs[g][i])));
394 bNonSelfExcl = FALSE0;
395 for (j = 0; j < isize[g+1]; j++)
396 {
397 jx = index[g+1][j];
398 if (!bExcl[jx])
399 {
400 pairs[g][i][k++] = jx;
401 }
402 else if (ix != jx)
403 {
404 /* Check if we have exclusions other than self exclusions */
405 bNonSelfExcl = TRUE1;
406 }
407 }
408 if (bNonSelfExcl)
409 {
410 npairs[g][i] = k;
411 srenew(pairs[g][i], npairs[g][i])(pairs[g][i]) = save_realloc("pairs[g][i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 411, (pairs[g][i]), (npairs[g][i]), sizeof(*(pairs[g][i])));
412 }
413 else
414 {
415 /* Save a LOT of memory and some cpu cycles */
416 npairs[g][i] = -1;
417 sfree(pairs[g][i])save_free("pairs[g][i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 417, (pairs[g][i]));
418 }
419 }
420 else
421 {
422 npairs[g][i] = -1;
423 }
424 }
425 }
426 sfree(bExcl)save_free("bExcl", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 426, (bExcl));
427
428 snew(x_i1, max_i)(x_i1) = save_calloc("x_i1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 428, (max_i), sizeof(*(x_i1)));
429 nframes = 0;
430 invvol_sum = 0;
431 if (bPBC && (NULL((void*)0) != top))
432 {
433 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
434 }
435
436 do
437 {
438 /* Must init pbc every step because of pressure coupling */
439 copy_mat(box, box_pbc);
440 if (bPBC)
441 {
442 if (top != NULL((void*)0))
443 {
444 gmx_rmpbc(gpbc, natoms, box, x);
445 }
446 if (bXY)
447 {
448 check_box_c(box);
449 clear_rvec(box_pbc[ZZ2]);
450 }
451 set_pbc(&pbc, ePBCrdf, box_pbc);
452
453 if (bXY)
454 {
455 /* Set z-size to 1 so we get the surface iso the volume */
456 box_pbc[ZZ2][ZZ2] = 1;
457 }
458 }
459 invvol = 1/det(box_pbc);
460 invvol_sum += invvol;
461
462 if (bCM)
463 {
464 /* Calculate center of mass of the whole group */
465 calc_comg(is[0], coi[0], index[0], TRUE1, atom, x, x0);
466 }
467 else if (!bClose && rdft[0][0] != 'a')
468 {
469 calc_comg(is[0], coi[0], index[0], rdft[0][6] == 'm', atom, x, x0);
470 }
471
472 for (g = 0; g < ng; g++)
473 {
474 if (rdft[0][0] == 'a')
475 {
476 /* Copy the indexed coordinates to a continuous array */
477 for (i = 0; i < isize[g+1]; i++)
478 {
479 copy_rvec(x[index[g+1][i]], x_i1[i]);
480 }
481 }
482 else
483 {
484 /* Calculate the COMs/COGs and store in x_i1 */
485 calc_comg(is[g+1], coi[g+1], index[g+1], rdft[0][6] == 'm', atom, x, x_i1);
486 }
487
488 for (i = 0; i < isize0; i++)
489 {
490 if (bClose)
491 {
492 /* Special loop, since we need to determine the minimum distance
493 * over all selected atoms in the reference molecule/residue.
494 */
495 if (rdft[0][0] == 'a')
496 {
497 isize_g = isize[g+1];
498 }
499 else
500 {
501 isize_g = is[g+1];
502 }
503 for (j = 0; j < isize_g; j++)
504 {
505 r2 = 1e30;
506 /* Loop over the selected atoms in the reference molecule */
507 for (ii = coi[0][i]; ii < coi[0][i+1]; ii++)
508 {
509 if (bPBC)
510 {
511 pbc_dx(&pbc, x[index[0][ii]], x_i1[j], dx);
512 }
513 else
514 {
515 rvec_sub(x[index[0][ii]], x_i1[j], dx);
516 }
517 if (bXY)
518 {
519 r2ii = dx[XX0]*dx[XX0] + dx[YY1]*dx[YY1];
520 }
521 else
522 {
523 r2ii = iprod(dx, dx);
524 }
525 if (r2ii < r2)
526 {
527 r2 = r2ii;
528 }
529 }
530 if (r2 > cut2 && r2 <= rmax2)
531 {
532 count[g][(int)(sqrt(r2)*invhbinw)]++;
533 }
534 }
535 }
536 else
537 {
538 /* Real rdf between points in space */
539 if (bCM || rdft[0][0] != 'a')
540 {
541 copy_rvec(x0[i], xi);
542 }
543 else
544 {
545 copy_rvec(x[index[0][i]], xi);
546 }
547 if (rdft[0][0] == 'a' && npairs[g][i] >= 0)
548 {
549 /* Expensive loop, because of indexing */
550 for (j = 0; j < npairs[g][i]; j++)
551 {
552 jx = pairs[g][i][j];
553 if (bPBC)
554 {
555 pbc_dx(&pbc, xi, x[jx], dx);
556 }
557 else
558 {
559 rvec_sub(xi, x[jx], dx);
560 }
561
562 if (bXY)
563 {
564 r2 = dx[XX0]*dx[XX0] + dx[YY1]*dx[YY1];
565 }
566 else
567 {
568 r2 = iprod(dx, dx);
569 }
570 if (r2 > cut2 && r2 <= rmax2)
571 {
572 count[g][(int)(sqrt(r2)*invhbinw)]++;
573 }
574 }
575 }
576 else
577 {
578 /* Cheaper loop, no exclusions */
579 if (rdft[0][0] == 'a')
580 {
581 isize_g = isize[g+1];
582 }
583 else
584 {
585 isize_g = is[g+1];
586 }
587 for (j = 0; j < isize_g; j++)
588 {
589 if (bPBC)
590 {
591 pbc_dx(&pbc, xi, x_i1[j], dx);
592 }
593 else
594 {
595 rvec_sub(xi, x_i1[j], dx);
596 }
597 if (bXY)
598 {
599 r2 = dx[XX0]*dx[XX0] + dx[YY1]*dx[YY1];
600 }
601 else
602 {
603 r2 = iprod(dx, dx);
604 }
605 if (r2 > cut2 && r2 <= rmax2)
606 {
607 count[g][(int)(sqrt(r2)*invhbinw)]++;
608 }
609 }
610 }
611 }
612 }
613 }
614 nframes++;
615 }
616 while (read_next_x(oenv, status, &t, x, box));
617 fprintf(stderrstderr, "\n");
618
619 if (bPBC && (NULL((void*)0) != top))
620 {
621 gmx_rmpbc_done(gpbc);
622 }
623
624 close_trj(status);
625
626 sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 626, (x));
627
628 /* Average volume */
629 invvol = invvol_sum/nframes;
630
631 /* Calculate volume of sphere segments or length of circle segments */
632 snew(inv_segvol, (nbin+1)/2)(inv_segvol) = save_calloc("inv_segvol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 632, ((nbin+1)/2), sizeof(*(inv_segvol)));
633 prev_spherevol = 0;
634 for (i = 0; (i < (nbin+1)/2); i++)
635 {
636 r = (i + 0.5)*binwidth;
637 if (bXY)
638 {
639 spherevol = M_PI3.14159265358979323846*r*r;
640 }
641 else
642 {
643 spherevol = (4.0/3.0)*M_PI3.14159265358979323846*r*r*r;
644 }
645 segvol = spherevol-prev_spherevol;
646 inv_segvol[i] = 1.0/segvol;
647 prev_spherevol = spherevol;
648 }
649
650 snew(rdf, ng)(rdf) = save_calloc("rdf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 650, (ng), sizeof(*(rdf)));
651 for (g = 0; g < ng; g++)
652 {
653 /* We have to normalize by dividing by the number of frames */
654 if (rdft[0][0] == 'a')
655 {
656 normfac = 1.0/(nframes*invvol*isize0*isize[g+1]);
657 }
658 else
659 {
660 normfac = 1.0/(nframes*invvol*isize0*is[g+1]);
661 }
662
663 /* Do the normalization */
664 nrdf = max((nbin+1)/2, 1+2*fade/binwidth)((((nbin+1)/2) > (1+2*fade/binwidth)) ? ((nbin+1)/2) : (1+
2*fade/binwidth) );
665 snew(rdf[g], nrdf)(rdf[g]) = save_calloc("rdf[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 665, (nrdf), sizeof(*(rdf[g])));
666 for (i = 0; i < (nbin+1)/2; i++)
667 {
668 r = i*binwidth;
669 if (i == 0)
670 {
671 j = count[g][0];
672 }
673 else
674 {
675 j = count[g][i*2-1] + count[g][i*2];
676 }
677 if ((fade > 0) && (r >= fade))
678 {
679 rdf[g][i] = 1 + (j*inv_segvol[i]*normfac-1)*exp(-16*sqr(r/fade-1));
680 }
681 else
682 {
683 if (bNormalize)
684 {
685 rdf[g][i] = j*inv_segvol[i]*normfac;
686 }
687 else
688 {
689 rdf[g][i] = j/(binwidth*isize0*nframes);
690 }
691 }
692 }
693 for (; (i < nrdf); i++)
694 {
695 rdf[g][i] = 1.0;
696 }
697 }
698
699 if (rdft[0][0] == 'a')
700 {
701 sprintf(gtitle, "Radial distribution");
702 }
703 else
704 {
705 sprintf(gtitle, "Radial distribution of %s %s",
706 rdft[0][0] == 'm' ? "molecule" : "residue",
707 rdft[0][6] == 'm' ? "COM" : "COG");
708 }
709 fp = xvgropen(fnRDF, gtitle, "r", "", oenv);
710 if (bCM)
711 {
712 sprintf(refgt, " %s", "COM");
713 }
714 else if (bClose)
715 {
716 sprintf(refgt, " closest atom in %s.", close);
717 }
718 else
719 {
720 sprintf(refgt, "%s", "");
721 }
722 if (ng == 1)
723 {
724 if (output_env_get_print_xvgr_codes(oenv))
725 {
726 fprintf(fp, "@ subtitle \"%s%s - %s\"\n", grpname[0], refgt, grpname[1]);
727 }
728 }
729 else
730 {
731 if (output_env_get_print_xvgr_codes(oenv))
732 {
733 fprintf(fp, "@ subtitle \"reference %s%s\"\n", grpname[0], refgt);
734 }
735 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
736 }
737 for (i = 0; (i < nrdf); i++)
738 {
739 fprintf(fp, "%10g", i*binwidth);
740 for (g = 0; g < ng; g++)
741 {
742 fprintf(fp, " %10g", rdf[g][i]);
743 }
744 fprintf(fp, "\n");
745 }
746 gmx_ffclose(fp);
747
748 do_view(oenv, fnRDF, NULL((void*)0));
749
750 /* h(Q) function: fourier transform of rdf */
751 if (fnHQ)
752 {
753 int nhq = 401;
754 real *hq, *integrand, Q;
755
756 /* Get a better number density later! */
757 rho = isize[1]*invvol;
758 snew(hq, nhq)(hq) = save_calloc("hq", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 758, (nhq), sizeof(*(hq)));
759 snew(integrand, nrdf)(integrand) = save_calloc("integrand", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 759, (nrdf), sizeof(*(integrand)));
760 for (i = 0; (i < nhq); i++)
761 {
762 Q = i*0.5;
763 integrand[0] = 0;
764 for (j = 1; (j < nrdf); j++)
765 {
766 r = j*binwidth;
767 integrand[j] = (Q == 0) ? 1.0 : sin(Q*r)/(Q*r);
768 integrand[j] *= 4.0*M_PI3.14159265358979323846*rho*r*r*(rdf[0][j]-1.0);
769 }
770 hq[i] = print_and_integrate(debug, nrdf, binwidth, integrand, NULL((void*)0), 0);
771 }
772 fp = xvgropen(fnHQ, "h(Q)", "Q(/nm)", "h(Q)", oenv);
773 for (i = 0; (i < nhq); i++)
774 {
775 fprintf(fp, "%10g %10g\n", i*0.5, hq[i]);
776 }
777 gmx_ffclose(fp);
778 do_view(oenv, fnHQ, NULL((void*)0));
779 sfree(hq)save_free("hq", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 779, (hq));
780 sfree(integrand)save_free("integrand", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 780, (integrand));
781 }
782
783 if (fnCNRDF)
784 {
785 normfac = 1.0/(isize0*nframes);
786 fp = xvgropen(fnCNRDF, "Cumulative Number RDF", "r", "number", oenv);
787 if (ng == 1)
788 {
789 if (output_env_get_print_xvgr_codes(oenv))
790 {
791 fprintf(fp, "@ subtitle \"%s-%s\"\n", grpname[0], grpname[1]);
792 }
793 }
794 else
795 {
796 if (output_env_get_print_xvgr_codes(oenv))
797 {
798 fprintf(fp, "@ subtitle \"reference %s\"\n", grpname[0]);
799 }
800 xvgr_legend(fp, ng, (const char**)(grpname+1), oenv);
801 }
802 snew(sum, ng)(sum) = save_calloc("sum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 802, (ng), sizeof(*(sum)));
803 for (i = 0; (i <= nbin/2); i++)
804 {
805 fprintf(fp, "%10g", i*binwidth);
806 for (g = 0; g < ng; g++)
807 {
808 fprintf(fp, " %10g", (real)((double)sum[g]*normfac));
809 if (i*2+1 < nbin)
810 {
811 sum[g] += count[g][i*2] + count[g][i*2+1];
812 }
813 }
814 fprintf(fp, "\n");
815 }
816 gmx_ffclose(fp);
817 sfree(sum)save_free("sum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 817, (sum));
818
819 do_view(oenv, fnCNRDF, NULL((void*)0));
820 }
821
822 for (g = 0; g < ng; g++)
823 {
824 sfree(rdf[g])save_free("rdf[g]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 824, (rdf[g]));
825 }
826 sfree(rdf)save_free("rdf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 826, (rdf));
827}
828
829
830int gmx_rdf(int argc, char *argv[])
831{
832 const char *desc[] = {
833 "The structure of liquids can be studied by either neutron or X-ray",
834 "scattering. The most common way to describe liquid structure is by a",
835 "radial distribution function. However, this is not easy to obtain from",
836 "a scattering experiment.[PAR]",
837 "[THISMODULE] calculates radial distribution functions in different ways.",
838 "The normal method is around a (set of) particle(s), the other methods",
839 "are around the center of mass of a set of particles ([TT]-com[tt])",
840 "or to the closest particle in a set ([TT]-surf[tt]).",
841 "With all methods, the RDF can also be calculated around axes parallel",
842 "to the [IT]z[it]-axis with option [TT]-xy[tt].",
843 "With option [TT]-surf[tt] normalization can not be used.[PAR]",
844 "The option [TT]-rdf[tt] sets the type of RDF to be computed.",
845 "Default is for atoms or particles, but one can also select center",
846 "of mass or geometry of molecules or residues. In all cases, only",
847 "the atoms in the index groups are taken into account.",
848 "For molecules and/or the center of mass option, a run input file",
849 "is required.",
850 "Weighting other than COM or COG can currently only be achieved",
851 "by providing a run input file with different masses.",
852 "Options [TT]-com[tt] and [TT]-surf[tt] also work in conjunction",
853 "with [TT]-rdf[tt].[PAR]",
854 "If a run input file is supplied ([TT]-s[tt]) and [TT]-rdf[tt] is set",
855 "to [TT]atom[tt], exclusions defined",
856 "in that file are taken into account when calculating the RDF.",
857 "The option [TT]-cut[tt] is meant as an alternative way to avoid",
858 "intramolecular peaks in the RDF plot.",
859 "It is however better to supply a run input file with a higher number of",
860 "exclusions. For e.g. benzene a topology, setting nrexcl to 5",
861 "would eliminate all intramolecular contributions to the RDF.",
862 "Note that all atoms in the selected groups are used, also the ones",
863 "that don't have Lennard-Jones interactions.[PAR]",
864 "Option [TT]-cn[tt] produces the cumulative number RDF,",
865 "i.e. the average number of particles within a distance r.[PAR]"
866 };
867 static gmx_bool bCM = FALSE0, bXY = FALSE0, bPBC = TRUE1, bNormalize = TRUE1;
868 static real cutoff = 0, binwidth = 0.002, fade = 0.0;
869 static int ngroups = 1;
870
871 static const char *closet[] = { NULL((void*)0), "no", "mol", "res", NULL((void*)0) };
872 static const char *rdft[] = { NULL((void*)0), "atom", "mol_com", "mol_cog", "res_com", "res_cog", NULL((void*)0) };
873
874 t_pargs pa[] = {
875 { "-bin", FALSE0, etREAL, {&binwidth},
876 "Binwidth (nm)" },
877 { "-com", FALSE0, etBOOL, {&bCM},
878 "RDF with respect to the center of mass of first group" },
879 { "-surf", FALSE0, etENUM, {closet},
880 "RDF with respect to the surface of the first group" },
881 { "-rdf", FALSE0, etENUM, {rdft},
882 "RDF type" },
883 { "-pbc", FALSE0, etBOOL, {&bPBC},
884 "Use periodic boundary conditions for computing distances. Without PBC the maximum range will be three times the largest box edge." },
885 { "-norm", FALSE0, etBOOL, {&bNormalize},
886 "Normalize for volume and density" },
887 { "-xy", FALSE0, etBOOL, {&bXY},
888 "Use only the x and y components of the distance" },
889 { "-cut", FALSE0, etREAL, {&cutoff},
890 "Shortest distance (nm) to be considered"},
891 { "-ng", FALSE0, etINT, {&ngroups},
892 "Number of secondary groups to compute RDFs around a central group" },
893 { "-fade", FALSE0, etREAL, {&fade},
894 "From this distance onwards the RDF is tranformed by g'(r) = 1 + [g(r)-1] exp(-(r/fade-1)^2 to make it go to 1 smoothly. If fade is 0.0 nothing is done." }
895 };
896#define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])))
897 const char *fnTPS, *fnNDX;
898 output_env_t oenv;
899
900 t_filenm fnm[] = {
901 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
902 { efTPS, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
903 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
904 { efXVG, "-o", "rdf", ffWRITE1<<2 },
905 { efXVG, "-cn", "rdf_cn", ffOPTWR(1<<2| 1<<3) },
906 { efXVG, "-hq", "hq", ffOPTWR(1<<2| 1<<3) },
907 };
908#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
909 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
→
910 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv))
911 {
912 return 0;
913 }
914
915 if (bCM || closet[0][0] != 'n' || rdft[0][0] == 'm' || rdft[0][6] == 'm')
2
←
Array access results in a null pointer dereference
916 {
917 fnTPS = ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
918 }
919 else
920 {
921 fnTPS = ftp2fn_null(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
922 }
923 fnNDX = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
924
925 if (!fnTPS && !fnNDX)
926 {
927 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 927, "Neither index file nor topology file specified\n"
928 "Nothing to do!");
929 }
930
931 if (closet[0][0] != 'n')
932 {
933 if (bCM)
934 {
935 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rdf.c"
, 935, "Can not have both -com and -surf");
936 }
937 if (bNormalize)
938 {
939 fprintf(stderrstderr, "Turning of normalization because of option -surf\n");
940 bNormalize = FALSE0;
941 }
942 }
943
944 do_rdf(fnNDX, fnTPS, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
945 opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn_null("-cn", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
946 opt2fn_null("-hq", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
947 bCM, closet[0], rdft, bXY, bPBC, bNormalize, cutoff, binwidth, fade, ngroups,
948 oenv);
949
950 return 0;
951}