File: | gromacs/gmxana/gmx_densmap.c |
Location: | line 239, column 5 |
Description: | Value stored to 'natoms' is never read |
1 | /* |
2 | * This file is part of the GROMACS molecular simulation package. |
3 | * |
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. |
5 | * Copyright (c) 2001-2004, The GROMACS development team. |
6 | * Copyright (c) 2013,2014, by the GROMACS development team, led by |
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, |
8 | * and including many others, as listed in the AUTHORS file in the |
9 | * top-level source directory and at http://www.gromacs.org. |
10 | * |
11 | * GROMACS is free software; you can redistribute it and/or |
12 | * modify it under the terms of the GNU Lesser General Public License |
13 | * as published by the Free Software Foundation; either version 2.1 |
14 | * of the License, or (at your option) any later version. |
15 | * |
16 | * GROMACS is distributed in the hope that it will be useful, |
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
19 | * Lesser General Public License for more details. |
20 | * |
21 | * You should have received a copy of the GNU Lesser General Public |
22 | * License along with GROMACS; if not, see |
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, |
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. |
25 | * |
26 | * If you want to redistribute modifications to GROMACS, please |
27 | * consider that scientific software is very special. Version |
28 | * control is crucial - bugs must be traceable. We will be happy to |
29 | * consider code for inclusion in the official distribution, but |
30 | * derived work must not be called official GROMACS. Details are found |
31 | * in the README & COPYING files - if they are missing, get the |
32 | * official version at http://www.gromacs.org. |
33 | * |
34 | * To help us fund GROMACS development, we humbly ask that you cite |
35 | * the research papers on the package. Check out http://www.gromacs.org. |
36 | */ |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <math.h> |
42 | #include <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 | |
66 | int 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), |
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) |
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) |
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)) |
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) |
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) |
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]) |
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); |
Value stored to 'natoms' is never read | |
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 | } |