Bug Summary

File:gromacs/gmxana/gmx_densmap.c
Location:line 232, 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
41#include <math.h>
42#include <string.h>
43
44#include "gromacs/commandline/pargs.h"
45#include "typedefs.h"
46#include "gromacs/utility/smalloc.h"
47#include "macros.h"
48#include "gromacs/math/vec.h"
49#include "pbc.h"
50#include "gromacs/utility/futil.h"
51#include "index.h"
52#include "mshift.h"
53#include "princ.h"
54#include "rmpbc.h"
55#include "txtdump.h"
56#include "gromacs/fileio/tpxio.h"
57#include "gromacs/fileio/trxio.h"
58#include "gstat.h"
59#include "gromacs/fileio/matio.h"
60#include "pbc.h"
61#include "viewit.h"
62#include "gmx_ana.h"
63
64#include "gromacs/utility/fatalerror.h"
65
66int gmx_densmap(int argc, char *argv[])
67{
68 const char *desc[] = {
69 "[THISMODULE] computes 2D number-density maps.",
70 "It can make planar and axial-radial density maps.",
71 "The output [TT].xpm[tt] file can be visualized with for instance xv",
72 "and can be converted to postscript with [TT]xpm2ps[tt].",
73 "Optionally, output can be in text form to a [TT].dat[tt] file with [TT]-od[tt], instead of the usual [TT].xpm[tt] file with [TT]-o[tt].",
74 "[PAR]",
75 "The default analysis is a 2-D number-density map for a selected",
76 "group of atoms in the x-y plane.",
77 "The averaging direction can be changed with the option [TT]-aver[tt].",
78 "When [TT]-xmin[tt] and/or [TT]-xmax[tt] are set only atoms that are",
79 "within the limit(s) in the averaging direction are taken into account.",
80 "The grid spacing is set with the option [TT]-bin[tt].",
81 "When [TT]-n1[tt] or [TT]-n2[tt] is non-zero, the grid",
82 "size is set by this option.",
83 "Box size fluctuations are properly taken into account.",
84 "[PAR]",
85 "When options [TT]-amax[tt] and [TT]-rmax[tt] are set, an axial-radial",
86 "number-density map is made. Three groups should be supplied, the centers",
87 "of mass of the first two groups define the axis, the third defines the",
88 "analysis group. The axial direction goes from -amax to +amax, where",
89 "the center is defined as the midpoint between the centers of mass and",
90 "the positive direction goes from the first to the second center of mass.",
91 "The radial direction goes from 0 to rmax or from -rmax to +rmax",
92 "when the [TT]-mirror[tt] option has been set.",
93 "[PAR]",
94 "The normalization of the output is set with the [TT]-unit[tt] option.",
95 "The default produces a true number density. Unit [TT]nm-2[tt] leaves out",
96 "the normalization for the averaging or the angular direction.",
97 "Option [TT]count[tt] produces the count for each grid cell.",
98 "When you do not want the scale in the output to go",
99 "from zero to the maximum density, you can set the maximum",
100 "with the option [TT]-dmax[tt]."
101 };
102 static int n1 = 0, n2 = 0;
103 static real xmin = -1, xmax = -1, bin = 0.02, dmin = 0, dmax = 0, amax = 0, rmax = 0;
104 static gmx_bool bMirror = FALSE0, bSums = FALSE0;
105 static const char *eaver[] = { NULL((void*)0), "z", "y", "x", NULL((void*)0) };
106 static const char *eunit[] = { NULL((void*)0), "nm-3", "nm-2", "count", NULL((void*)0) };
107
108 t_pargs pa[] = {
109 { "-bin", FALSE0, etREAL, {&bin},
110 "Grid size (nm)" },
111 { "-aver", FALSE0, etENUM, {eaver},
112 "The direction to average over" },
113 { "-xmin", FALSE0, etREAL, {&xmin},
114 "Minimum coordinate for averaging" },
115 { "-xmax", FALSE0, etREAL, {&xmax},
116 "Maximum coordinate for averaging" },
117 { "-n1", FALSE0, etINT, {&n1},
118 "Number of grid cells in the first direction" },
119 { "-n2", FALSE0, etINT, {&n2},
120 "Number of grid cells in the second direction" },
121 { "-amax", FALSE0, etREAL, {&amax},
122 "Maximum axial distance from the center"},
123 { "-rmax", FALSE0, etREAL, {&rmax},
124 "Maximum radial distance" },
125 { "-mirror", FALSE0, etBOOL, {&bMirror},
126 "Add the mirror image below the axial axis" },
127 { "-sums", FALSE0, etBOOL, {&bSums},
128 "Print density sums (1D map) to stdout" },
129 { "-unit", FALSE0, etENUM, {eunit},
130 "Unit for the output" },
131 { "-dmin", FALSE0, etREAL, {&dmin},
132 "Minimum density in output"},
133 { "-dmax", FALSE0, etREAL, {&dmax},
134 "Maximum density in output (0 means calculate it)"},
135 };
136 gmx_bool bXmin, bXmax, bRadial;
137 FILE *fp;
138 t_trxstatus *status;
139 t_topology top;
140 int ePBC = -1;
141 rvec *x, xcom[2], direction, center, dx;
142 matrix box;
143 real t, m, mtot;
144 t_pbc pbc;
145 int cav = 0, c1 = 0, c2 = 0, natoms;
146 char **grpname, title[256], buf[STRLEN4096];
147 const char *unit;
148 int i, j, k, l, ngrps, anagrp, *gnx = NULL((void*)0), nindex, nradial = 0, nfr, nmpower;
149 atom_id **ind = NULL((void*)0), *index;
150 real **grid, maxgrid, m1, m2, box1, box2, *tickx, *tickz, invcellvol;
151 real invspa = 0, invspz = 0, axial, r, vol_old, vol, rowsum;
152 int nlev = 51;
153 t_rgb rlo = {1, 1, 1}, rhi = {0, 0, 0};
154 output_env_t oenv;
155 const char *label[] = { "x (nm)", "y (nm)", "z (nm)" };
156 t_filenm fnm[] = {
157 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
158 { efTPS, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
159 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
160 { efDAT, "-od", "densmap", ffOPTWR(1<<2| 1<<3) },
161 { efXPM, "-o", "densmap", ffWRITE1<<2 }
162 };
163#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
164 int npargs;
165
166 npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])));
167
168 if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_CAN_VIEW(1<<5) | PCA_BE_NICE(1<<13),
1
Taking false branch
→
169 NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv))
170 {
171 return 0;
172 }
173
174 bXmin = opt2parg_bSet("-xmin", npargs, pa);
175 bXmax = opt2parg_bSet("-xmax", npargs, pa);
176 bRadial = (amax > 0 || rmax > 0);
177 if (bRadial)
2
←
Assuming 'bRadial' is 0
→
3
←
Taking false branch
→
178 {
179 if (amax <= 0 || rmax <= 0)
180 {
181 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 181, "Both amax and rmax should be larger than zero");
182 }
183 }
184
185 if (strcmp(eunit[0], "nm-3")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(eunit[0]) && __builtin_constant_p ("nm-3") &&
(__s1_len = strlen (eunit[0]), __s2_len = strlen ("nm-3"), (
!((size_t)(const void *)((eunit[0]) + 1) - (size_t)(const void
*)(eunit[0]) == 1) || __s1_len >= 4) && (!((size_t
)(const void *)(("nm-3") + 1) - (size_t)(const void *)("nm-3"
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (eunit[0], "nm-3"
) : (__builtin_constant_p (eunit[0]) && ((size_t)(const
void *)((eunit[0]) + 1) - (size_t)(const void *)(eunit[0]) ==
1) && (__s1_len = strlen (eunit[0]), __s1_len < 4
) ? (__builtin_constant_p ("nm-3") && ((size_t)(const
void *)(("nm-3") + 1) - (size_t)(const void *)("nm-3") == 1)
? __builtin_strcmp (eunit[0], "nm-3") : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
("nm-3"); int __result = (((const unsigned char *) (const char
*) (eunit[0]))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (eunit[0]))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (eunit[0]))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (eunit[0]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("nm-3") && ((size_t)(const void *)(("nm-3") + 1) - (
size_t)(const void *)("nm-3") == 1) && (__s2_len = strlen
("nm-3"), __s2_len < 4) ? (__builtin_constant_p (eunit[0]
) && ((size_t)(const void *)((eunit[0]) + 1) - (size_t
)(const void *)(eunit[0]) == 1) ? __builtin_strcmp (eunit[0],
"nm-3") : (- (__extension__ ({ const unsigned char *__s2 = (
const unsigned char *) (const char *) (eunit[0]); int __result
= (((const unsigned char *) (const char *) ("nm-3"))[0] - __s2
[0]); if (__s2_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("nm-3"))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("nm-3"))[2] - __s2
[2]); if (__s2_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) ("nm-3"))[3] - __s2
[3]); } } __result; })))) : __builtin_strcmp (eunit[0], "nm-3"
)))); }) == 0)
4
←
Taking true branch
→
186 {
187 nmpower = -3;
188 unit = "(nm^-3)";
189 }
190 else if (strcmp(eunit[0], "nm-2")__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(eunit[0]) && __builtin_constant_p ("nm-2") &&
(__s1_len = strlen (eunit[0]), __s2_len = strlen ("nm-2"), (
!((size_t)(const void *)((eunit[0]) + 1) - (size_t)(const void
*)(eunit[0]) == 1) || __s1_len >= 4) && (!((size_t
)(const void *)(("nm-2") + 1) - (size_t)(const void *)("nm-2"
) == 1) || __s2_len >= 4)) ? __builtin_strcmp (eunit[0], "nm-2"
) : (__builtin_constant_p (eunit[0]) && ((size_t)(const
void *)((eunit[0]) + 1) - (size_t)(const void *)(eunit[0]) ==
1) && (__s1_len = strlen (eunit[0]), __s1_len < 4
) ? (__builtin_constant_p ("nm-2") && ((size_t)(const
void *)(("nm-2") + 1) - (size_t)(const void *)("nm-2") == 1)
? __builtin_strcmp (eunit[0], "nm-2") : (__extension__ ({ const
unsigned char *__s2 = (const unsigned char *) (const char *)
("nm-2"); int __result = (((const unsigned char *) (const char
*) (eunit[0]))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (eunit[0]))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (eunit[0]))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (eunit[0]))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
("nm-2") && ((size_t)(const void *)(("nm-2") + 1) - (
size_t)(const void *)("nm-2") == 1) && (__s2_len = strlen
("nm-2"), __s2_len < 4) ? (__builtin_constant_p (eunit[0]
) && ((size_t)(const void *)((eunit[0]) + 1) - (size_t
)(const void *)(eunit[0]) == 1) ? __builtin_strcmp (eunit[0],
"nm-2") : (- (__extension__ ({ const unsigned char *__s2 = (
const unsigned char *) (const char *) (eunit[0]); int __result
= (((const unsigned char *) (const char *) ("nm-2"))[0] - __s2
[0]); if (__s2_len > 0 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("nm-2"))[1] - __s2
[1]); if (__s2_len > 1 && __result == 0) { __result
= (((const unsigned char *) (const char *) ("nm-2"))[2] - __s2
[2]); if (__s2_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) ("nm-2"))[3] - __s2
[3]); } } __result; })))) : __builtin_strcmp (eunit[0], "nm-2"
)))); }) == 0)
191 {
192 nmpower = -2;
193 unit = "(nm^-2)";
194 }
195 else
196 {
197 nmpower = 0;
198 unit = "count";
199 }
200
201 if (ftp2bSet(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || !ftp2bSet(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
5
←
Taking false branch
→
202 {
203 read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &x, NULL((void*)0), box,
204 bRadial);
205 }
206 if (!bRadial)
6
←
Taking true branch
→
207 {
208 ngrps = 1;
209 fprintf(stderrstderr, "\nSelect an analysis group\n");
210 }
211 else
212 {
213 ngrps = 3;
214 fprintf(stderrstderr,
215 "\nSelect two groups to define the axis and an analysis group\n");
216 }
217 snew(gnx, ngrps)(gnx) = save_calloc("gnx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 217, (ngrps), sizeof(*(gnx)));
218 snew(grpname, ngrps)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 218, (ngrps), sizeof(*(grpname)));
219 snew(ind, ngrps)(ind) = save_calloc("ind", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 219, (ngrps), sizeof(*(ind)));
220 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), ngrps, gnx, ind, grpname);
221 anagrp = ngrps - 1;
222 nindex = gnx[anagrp];
223 index = ind[anagrp];
224 if (bRadial)
7
←
Taking false branch
→
225 {
226 if ((gnx[0] > 1 || gnx[1] > 1) && !ftp2bSet(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
227 {
228 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 228, "No run input file was supplied (option -s), this is required for the center of mass calculation");
229 }
230 }
231
232 switch (eaver[0][0])
8
←
Array access results in a null pointer dereference
233 {
234 case 'x': cav = XX0; c1 = YY1; c2 = ZZ2; break;
235 case 'y': cav = YY1; c1 = XX0; c2 = ZZ2; break;
236 case 'z': cav = ZZ2; c1 = XX0; c2 = YY1; break;
237 }
238
239 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box);
240
241 if (!bRadial)
242 {
243 if (n1 == 0)
244 {
245 n1 = (int)(box[c1][c1]/bin + 0.5);
246 }
247 if (n2 == 0)
248 {
249 n2 = (int)(box[c2][c2]/bin + 0.5);
250 }
251 }
252 else
253 {
254 n1 = (int)(2*amax/bin + 0.5);
255 nradial = (int)(rmax/bin + 0.5);
256 invspa = n1/(2*amax);
257 invspz = nradial/rmax;
258 if (bMirror)
259 {
260 n2 = 2*nradial;
261 }
262 else
263 {
264 n2 = nradial;
265 }
266 }
267
268 snew(grid, n1)(grid) = save_calloc("grid", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 268, (n1), sizeof(*(grid)));
269 for (i = 0; i < n1; i++)
270 {
271 snew(grid[i], n2)(grid[i]) = save_calloc("grid[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 271, (n2), sizeof(*(grid[i])));
272 }
273
274 box1 = 0;
275 box2 = 0;
276 nfr = 0;
277 do
278 {
279 if (!bRadial)
280 {
281 box1 += box[c1][c1];
282 box2 += box[c2][c2];
283 invcellvol = n1*n2;
284 if (nmpower == -3)
285 {
286 invcellvol /= det(box);
287 }
288 else if (nmpower == -2)
289 {
290 invcellvol /= box[c1][c1]*box[c2][c2];
291 }
292 for (i = 0; i < nindex; i++)
293 {
294 j = index[i];
295 if ((!bXmin || x[j][cav] >= xmin) &&
296 (!bXmax || x[j][cav] <= xmax))
297 {
298 m1 = x[j][c1]/box[c1][c1];
299 if (m1 >= 1)
300 {
301 m1 -= 1;
302 }
303 if (m1 < 0)
304 {
305 m1 += 1;
306 }
307 m2 = x[j][c2]/box[c2][c2];
308 if (m2 >= 1)
309 {
310 m2 -= 1;
311 }
312 if (m2 < 0)
313 {
314 m2 += 1;
315 }
316 grid[(int)(m1*n1)][(int)(m2*n2)] += invcellvol;
317 }
318 }
319 }
320 else
321 {
322 set_pbc(&pbc, ePBC, box);
323 for (i = 0; i < 2; i++)
324 {
325 if (gnx[i] == 1)
326 {
327 /* One atom, just copy the coordinates */
328 copy_rvec(x[ind[i][0]], xcom[i]);
329 }
330 else
331 {
332 /* Calculate the center of mass */
333 clear_rvec(xcom[i]);
334 mtot = 0;
335 for (j = 0; j < gnx[i]; j++)
336 {
337 k = ind[i][j];
338 m = top.atoms.atom[k].m;
339 for (l = 0; l < DIM3; l++)
340 {
341 xcom[i][l] += m*x[k][l];
342 }
343 mtot += m;
344 }
345 svmul(1/mtot, xcom[i], xcom[i]);
346 }
347 }
348 pbc_dx(&pbc, xcom[1], xcom[0], direction);
349 for (i = 0; i < DIM3; i++)
350 {
351 center[i] = xcom[0][i] + 0.5*direction[i];
352 }
353 unitv(direction, direction);
354 for (i = 0; i < nindex; i++)
355 {
356 j = index[i];
357 pbc_dx(&pbc, x[j], center, dx);
358 axial = iprod(dx, direction);
359 r = sqrt(norm2(dx) - axial*axial);
360 if (axial >= -amax && axial < amax && r < rmax)
361 {
362 if (bMirror)
363 {
364 r += rmax;
365 }
366 grid[(int)((axial + amax)*invspa)][(int)(r*invspz)] += 1;
367 }
368 }
369 }
370 nfr++;
371 }
372 while (read_next_x(oenv, status, &t, x, box));
373 close_trj(status);
374
375 /* normalize gridpoints */
376 maxgrid = 0;
377 if (!bRadial)
378 {
379 for (i = 0; i < n1; i++)
380 {
381 for (j = 0; j < n2; j++)
382 {
383 grid[i][j] /= nfr;
384 if (grid[i][j] > maxgrid)
385 {
386 maxgrid = grid[i][j];
387 }
388 }
389 }
390 }
391 else
392 {
393 for (i = 0; i < n1; i++)
394 {
395 vol_old = 0;
396 for (j = 0; j < nradial; j++)
397 {
398 switch (nmpower)
399 {
400 case -3: vol = M_PI3.14159265358979323846*(j+1)*(j+1)/(invspz*invspz*invspa); break;
401 case -2: vol = (j+1)/(invspz*invspa); break;
402 default: vol = j+1; break;
403 }
404 if (bMirror)
405 {
406 k = j + nradial;
407 }
408 else
409 {
410 k = j;
411 }
412 grid[i][k] /= nfr*(vol - vol_old);
413 if (bMirror)
414 {
415 grid[i][nradial-1-j] = grid[i][k];
416 }
417 vol_old = vol;
418 if (grid[i][k] > maxgrid)
419 {
420 maxgrid = grid[i][k];
421 }
422 }
423 }
424 }
425 fprintf(stdoutstdout, "\n The maximum density is %f %s\n", maxgrid, unit);
426 if (dmax > 0)
427 {
428 maxgrid = dmax;
429 }
430
431 snew(tickx, n1+1)(tickx) = save_calloc("tickx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 431, (n1+1), sizeof(*(tickx)));
432 snew(tickz, n2+1)(tickz) = save_calloc("tickz", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_densmap.c"
, 432, (n2+1), sizeof(*(tickz)));
433 if (!bRadial)
434 {
435 /* normalize box-axes */
436 box1 /= nfr;
437 box2 /= nfr;
438 for (i = 0; i <= n1; i++)
439 {
440 tickx[i] = i*box1/n1;
441 }
442 for (i = 0; i <= n2; i++)
443 {
444 tickz[i] = i*box2/n2;
445 }
446 }
447 else
448 {
449 for (i = 0; i <= n1; i++)
450 {
451 tickx[i] = i/invspa - amax;
452 }
453 if (bMirror)
454 {
455 for (i = 0; i <= n2; i++)
456 {
457 tickz[i] = i/invspz - rmax;
458 }
459 }
460 else
461 {
462 for (i = 0; i <= n2; i++)
463 {
464 tickz[i] = i/invspz;
465 }
466 }
467 }
468
469 if (bSums)
470 {
471 for (i = 0; i < n1; ++i)
472 {
473 fprintf(stdoutstdout, "Density sums:\n");
474 rowsum = 0;
475 for (j = 0; j < n2; ++j)
476 {
477 rowsum += grid[i][j];
478 }
479 fprintf(stdoutstdout, "%g\t", rowsum);
480 }
481 fprintf(stdoutstdout, "\n");
482 }
483
484 sprintf(buf, "%s number density", grpname[anagrp]);
485 if (!bRadial && (bXmin || bXmax))
486 {
487 if (!bXmax)
488 {
489 sprintf(buf+strlen(buf), ", %c > %g nm", eaver[0][0], xmin);
490 }
491 else if (!bXmin)
492 {
493 sprintf(buf+strlen(buf), ", %c < %g nm", eaver[0][0], xmax);
494 }
495 else
496 {
497 sprintf(buf+strlen(buf), ", %c: %g - %g nm", eaver[0][0], xmin, xmax);
498 }
499 }
500 if (ftp2bSet(efDAT, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
501 {
502 fp = gmx_ffopen(ftp2fn(efDAT, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "w");
503 /*optional text form output: first row is tickz; first col is tickx */
504 fprintf(fp, "0\t");
505 for (j = 0; j < n2; ++j)
506 {
507 fprintf(fp, "%g\t", tickz[j]);
508 }
509 fprintf(fp, "\n");
510
511 for (i = 0; i < n1; ++i)
512 {
513 fprintf(fp, "%g\t", tickx[i]);
514 for (j = 0; j < n2; ++j)
515 {
516 fprintf(fp, "%g\t", grid[i][j]);
517 }
518 fprintf(fp, "\n");
519 }
520 gmx_ffclose(fp);
521 }
522 else
523 {
524 fp = gmx_ffopen(ftp2fn(efXPM, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "w");
525 write_xpm(fp, MAT_SPATIAL_X(1<<0) | MAT_SPATIAL_Y(1<<1), buf, unit,
526 bRadial ? "axial (nm)" : label[c1], bRadial ? "r (nm)" : label[c2],
527 n1, n2, tickx, tickz, grid, dmin, maxgrid, rlo, rhi, &nlev);
528 gmx_ffclose(fp);
529 }
530
531 do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
532
533 return 0;
534}