File: | gromacs/gmxana/gmx_density.c |
Location: | line 525, column 9 |
Description: | Array access results in a null pointer dereference |
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 <ctype.h> | |||
42 | #include <math.h> | |||
43 | #include <stdlib.h> | |||
44 | #include <string.h> | |||
45 | ||||
46 | #include "gromacs/utility/cstringutil.h" | |||
47 | #include "typedefs.h" | |||
48 | #include "gromacs/utility/smalloc.h" | |||
49 | #include "macros.h" | |||
50 | #include "gstat.h" | |||
51 | #include "gromacs/math/vec.h" | |||
52 | #include "gromacs/fileio/xvgr.h" | |||
53 | #include "viewit.h" | |||
54 | #include "pbc.h" | |||
55 | #include "gromacs/utility/futil.h" | |||
56 | #include "gromacs/commandline/pargs.h" | |||
57 | #include "index.h" | |||
58 | #include "gromacs/fileio/tpxio.h" | |||
59 | #include "gromacs/fileio/trxio.h" | |||
60 | #include "physics.h" | |||
61 | #include "gmx_ana.h" | |||
62 | #include "macros.h" | |||
63 | ||||
64 | #include "gromacs/utility/fatalerror.h" | |||
65 | ||||
66 | typedef struct { | |||
67 | char *atomname; | |||
68 | int nr_el; | |||
69 | } t_electron; | |||
70 | ||||
71 | /****************************************************************************/ | |||
72 | /* This program calculates the partial density across the box. */ | |||
73 | /* Peter Tieleman, Mei 1995 */ | |||
74 | /****************************************************************************/ | |||
75 | ||||
76 | /* used for sorting the list */ | |||
77 | int compare(void *a, void *b) | |||
78 | { | |||
79 | t_electron *tmp1, *tmp2; | |||
80 | tmp1 = (t_electron *)a; tmp2 = (t_electron *)b; | |||
81 | ||||
82 | return strcmp(tmp1->atomname, tmp2->atomname)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p (tmp1->atomname) && __builtin_constant_p (tmp2-> atomname) && (__s1_len = strlen (tmp1->atomname), __s2_len = strlen (tmp2->atomname), (!((size_t)(const void *)((tmp1 ->atomname) + 1) - (size_t)(const void *)(tmp1->atomname ) == 1) || __s1_len >= 4) && (!((size_t)(const void *)((tmp2->atomname) + 1) - (size_t)(const void *)(tmp2-> atomname) == 1) || __s2_len >= 4)) ? __builtin_strcmp (tmp1 ->atomname, tmp2->atomname) : (__builtin_constant_p (tmp1 ->atomname) && ((size_t)(const void *)((tmp1->atomname ) + 1) - (size_t)(const void *)(tmp1->atomname) == 1) && (__s1_len = strlen (tmp1->atomname), __s1_len < 4) ? ( __builtin_constant_p (tmp2->atomname) && ((size_t) (const void *)((tmp2->atomname) + 1) - (size_t)(const void *)(tmp2->atomname) == 1) ? __builtin_strcmp (tmp1->atomname , tmp2->atomname) : (__extension__ ({ const unsigned char * __s2 = (const unsigned char *) (const char *) (tmp2->atomname ); int __result = (((const unsigned char *) (const char *) (tmp1 ->atomname))[0] - __s2[0]); if (__s1_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (tmp1->atomname))[1] - __s2[1]); if (__s1_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) (tmp1->atomname))[2] - __s2[2]); if (__s1_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (tmp1->atomname))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (tmp2->atomname) && ((size_t)(const void *)((tmp2->atomname) + 1) - (size_t)(const void *)(tmp2->atomname) == 1) && ( __s2_len = strlen (tmp2->atomname), __s2_len < 4) ? (__builtin_constant_p (tmp1->atomname) && ((size_t)(const void *)((tmp1 ->atomname) + 1) - (size_t)(const void *)(tmp1->atomname ) == 1) ? __builtin_strcmp (tmp1->atomname, tmp2->atomname ) : (- (__extension__ ({ const unsigned char *__s2 = (const unsigned char *) (const char *) (tmp1->atomname); int __result = ( ((const unsigned char *) (const char *) (tmp2->atomname))[ 0] - __s2[0]); if (__s2_len > 0 && __result == 0) { __result = (((const unsigned char *) (const char *) (tmp2-> atomname))[1] - __s2[1]); if (__s2_len > 1 && __result == 0) { __result = (((const unsigned char *) (const char *) ( tmp2->atomname))[2] - __s2[2]); if (__s2_len > 2 && __result == 0) __result = (((const unsigned char *) (const char *) (tmp2->atomname))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (tmp1->atomname, tmp2->atomname)))); } ); | |||
83 | } | |||
84 | ||||
85 | int get_electrons(t_electron **eltab, const char *fn) | |||
86 | { | |||
87 | char buffer[256]; /* to read in a line */ | |||
88 | char tempname[80]; /* buffer to hold name */ | |||
89 | int tempnr; | |||
90 | ||||
91 | FILE *in; | |||
92 | int nr; /* number of atomstypes to read */ | |||
93 | int i; | |||
94 | ||||
95 | if (!(in = gmx_ffopen(fn, "r"))) | |||
96 | { | |||
97 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 97, "Couldn't open %s. Exiting.\n", fn); | |||
98 | } | |||
99 | ||||
100 | if (NULL((void*)0) == fgets(buffer, 255, in)) | |||
101 | { | |||
102 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 102, "Error reading from file %s", fn); | |||
103 | } | |||
104 | ||||
105 | if (sscanf(buffer, "%d", &nr) != 1) | |||
106 | { | |||
107 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 107, "Invalid number of atomtypes in datafile\n"); | |||
108 | } | |||
109 | ||||
110 | snew(*eltab, nr)(*eltab) = save_calloc("*eltab", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 110, (nr), sizeof(*(*eltab))); | |||
111 | ||||
112 | for (i = 0; i < nr; i++) | |||
113 | { | |||
114 | if (fgets(buffer, 255, in) == NULL((void*)0)) | |||
115 | { | |||
116 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 116, "reading datafile. Check your datafile.\n"); | |||
117 | } | |||
118 | if (sscanf(buffer, "%s = %d", tempname, &tempnr) != 2) | |||
119 | { | |||
120 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 120, "Invalid line in datafile at line %d\n", i+1); | |||
121 | } | |||
122 | (*eltab)[i].nr_el = tempnr; | |||
123 | (*eltab)[i].atomname = strdup(tempname)(__extension__ (__builtin_constant_p (tempname) && (( size_t)(const void *)((tempname) + 1) - (size_t)(const void * )(tempname) == 1) ? (((const char *) (tempname))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (tempname) + 1; char *__retval = (char *) malloc (__len ); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval , tempname, __len); __retval; })) : __strdup (tempname))); | |||
124 | } | |||
125 | gmx_ffclose(in); | |||
126 | ||||
127 | /* sort the list */ | |||
128 | fprintf(stderrstderr, "Sorting list..\n"); | |||
129 | qsort ((void*)*eltab, nr, sizeof(t_electron), | |||
130 | (int(*)(const void*, const void*))compare); | |||
131 | ||||
132 | return nr; | |||
133 | } | |||
134 | ||||
135 | void center_coords(t_atoms *atoms, matrix box, rvec x0[], int axis) | |||
136 | { | |||
137 | int i, m; | |||
138 | real tmass, mm; | |||
139 | rvec com, shift, box_center; | |||
140 | ||||
141 | tmass = 0; | |||
142 | clear_rvec(com); | |||
143 | for (i = 0; (i < atoms->nr); i++) | |||
144 | { | |||
145 | mm = atoms->atom[i].m; | |||
146 | tmass += mm; | |||
147 | for (m = 0; (m < DIM3); m++) | |||
148 | { | |||
149 | com[m] += mm*x0[i][m]; | |||
150 | } | |||
151 | } | |||
152 | for (m = 0; (m < DIM3); m++) | |||
153 | { | |||
154 | com[m] /= tmass; | |||
155 | } | |||
156 | calc_box_center(ecenterDEF, box, box_center); | |||
157 | rvec_sub(box_center, com, shift); | |||
158 | shift[axis] -= box_center[axis]; | |||
159 | ||||
160 | for (i = 0; (i < atoms->nr); i++) | |||
161 | { | |||
162 | rvec_dec(x0[i], shift); | |||
163 | } | |||
164 | } | |||
165 | ||||
166 | void calc_electron_density(const char *fn, atom_id **index, int gnx[], | |||
167 | double ***slDensity, int *nslices, t_topology *top, | |||
168 | int ePBC, | |||
169 | int axis, int nr_grps, real *slWidth, | |||
170 | t_electron eltab[], int nr, gmx_bool bCenter, | |||
171 | const output_env_t oenv) | |||
172 | { | |||
173 | rvec *x0; /* coordinates without pbc */ | |||
174 | matrix box; /* box (3x3) */ | |||
175 | double invvol; | |||
176 | int natoms; /* nr. atoms in trj */ | |||
177 | t_trxstatus *status; | |||
178 | int i, n, /* loop indices */ | |||
179 | nr_frames = 0, /* number of frames */ | |||
180 | slice; /* current slice */ | |||
181 | t_electron *found; /* found by bsearch */ | |||
182 | t_electron sought; /* thingie thought by bsearch */ | |||
183 | gmx_rmpbc_t gpbc = NULL((void*)0); | |||
184 | ||||
185 | real t, | |||
186 | z; | |||
187 | ||||
188 | if (axis < 0 || axis >= DIM3) | |||
189 | { | |||
190 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 190, "Invalid axes. Terminating\n"); | |||
191 | } | |||
192 | ||||
193 | if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0) | |||
194 | { | |||
195 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 195, "Could not read coordinates from statusfile\n"); | |||
196 | } | |||
197 | ||||
198 | if (!*nslices) | |||
199 | { | |||
200 | *nslices = (int)(box[axis][axis] * 10); /* default value */ | |||
201 | } | |||
202 | fprintf(stderrstderr, "\nDividing the box in %d slices\n", *nslices); | |||
203 | ||||
204 | snew(*slDensity, nr_grps)(*slDensity) = save_calloc("*slDensity", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 204, (nr_grps), sizeof(*(*slDensity))); | |||
205 | for (i = 0; i < nr_grps; i++) | |||
206 | { | |||
207 | snew((*slDensity)[i], *nslices)((*slDensity)[i]) = save_calloc("(*slDensity)[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 207, (*nslices), sizeof(*((*slDensity)[i]))); | |||
208 | } | |||
209 | ||||
210 | gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr); | |||
211 | /*********** Start processing trajectory ***********/ | |||
212 | do | |||
213 | { | |||
214 | gmx_rmpbc(gpbc, natoms, box, x0); | |||
215 | ||||
216 | if (bCenter) | |||
217 | { | |||
218 | center_coords(&top->atoms, box, x0, axis); | |||
219 | } | |||
220 | ||||
221 | *slWidth = box[axis][axis]/(*nslices); | |||
222 | invvol = *nslices/(box[XX0][XX0]*box[YY1][YY1]*box[ZZ2][ZZ2]); | |||
223 | ||||
224 | for (n = 0; n < nr_grps; n++) | |||
225 | { | |||
226 | for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */ | |||
227 | { | |||
228 | z = x0[index[n][i]][axis]; | |||
229 | while (z < 0) | |||
230 | { | |||
231 | z += box[axis][axis]; | |||
232 | } | |||
233 | while (z > box[axis][axis]) | |||
234 | { | |||
235 | z -= box[axis][axis]; | |||
236 | } | |||
237 | ||||
238 | /* determine which slice atom is in */ | |||
239 | slice = (z / (*slWidth)); | |||
240 | sought.nr_el = 0; | |||
241 | sought.atomname = strdup(*(top->atoms.atomname[index[n][i]]))(__extension__ (__builtin_constant_p (*(top->atoms.atomname [index[n][i]])) && ((size_t)(const void *)((*(top-> atoms.atomname[index[n][i]])) + 1) - (size_t)(const void *)(* (top->atoms.atomname[index[n][i]])) == 1) ? (((const char * ) (*(top->atoms.atomname[index[n][i]])))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (*(top->atoms.atomname[index[n][i]])) + 1; char *__retval = (char *) malloc (__len); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval, *(top->atoms.atomname[index[ n][i]]), __len); __retval; })) : __strdup (*(top->atoms.atomname [index[n][i]])))); | |||
242 | ||||
243 | /* now find the number of electrons. This is not efficient. */ | |||
244 | found = (t_electron *) | |||
245 | bsearch((const void *)&sought, | |||
246 | (const void *)eltab, nr, sizeof(t_electron), | |||
247 | (int(*)(const void*, const void*))compare); | |||
248 | ||||
249 | if (found == NULL((void*)0)) | |||
250 | { | |||
251 | fprintf(stderrstderr, "Couldn't find %s. Add it to the .dat file\n", | |||
252 | *(top->atoms.atomname[index[n][i]])); | |||
253 | } | |||
254 | else | |||
255 | { | |||
256 | (*slDensity)[n][slice] += (found->nr_el - | |||
257 | top->atoms.atom[index[n][i]].q)*invvol; | |||
258 | } | |||
259 | free(sought.atomname); | |||
260 | } | |||
261 | } | |||
262 | nr_frames++; | |||
263 | } | |||
264 | while (read_next_x(oenv, status, &t, x0, box)); | |||
265 | gmx_rmpbc_done(gpbc); | |||
266 | ||||
267 | /*********** done with status file **********/ | |||
268 | close_trj(status); | |||
269 | ||||
270 | /* slDensity now contains the total number of electrons per slice, summed | |||
271 | over all frames. Now divide by nr_frames and volume of slice | |||
272 | */ | |||
273 | ||||
274 | fprintf(stderrstderr, "\nRead %d frames from trajectory. Counting electrons\n", | |||
275 | nr_frames); | |||
276 | ||||
277 | for (n = 0; n < nr_grps; n++) | |||
278 | { | |||
279 | for (i = 0; i < *nslices; i++) | |||
280 | { | |||
281 | (*slDensity)[n][i] /= nr_frames; | |||
282 | } | |||
283 | } | |||
284 | ||||
285 | sfree(x0)save_free("x0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 285, (x0)); /* free memory used by coordinate array */ | |||
286 | } | |||
287 | ||||
288 | void calc_density(const char *fn, atom_id **index, int gnx[], | |||
289 | double ***slDensity, int *nslices, t_topology *top, int ePBC, | |||
290 | int axis, int nr_grps, real *slWidth, gmx_bool bCenter, | |||
291 | const output_env_t oenv) | |||
292 | { | |||
293 | rvec *x0; /* coordinates without pbc */ | |||
294 | matrix box; /* box (3x3) */ | |||
295 | double invvol; | |||
296 | int natoms; /* nr. atoms in trj */ | |||
297 | t_trxstatus *status; | |||
298 | int **slCount, /* nr. of atoms in one slice for a group */ | |||
299 | i, j, n, /* loop indices */ | |||
300 | teller = 0, | |||
301 | ax1 = 0, ax2 = 0, | |||
302 | nr_frames = 0, /* number of frames */ | |||
303 | slice; /* current slice */ | |||
304 | real t, | |||
305 | z; | |||
306 | char *buf; /* for tmp. keeping atomname */ | |||
307 | gmx_rmpbc_t gpbc = NULL((void*)0); | |||
308 | ||||
309 | if (axis < 0 || axis >= DIM3) | |||
310 | { | |||
311 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 311, "Invalid axes. Terminating\n"); | |||
312 | } | |||
313 | ||||
314 | if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0) | |||
315 | { | |||
316 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 316, "Could not read coordinates from statusfile\n"); | |||
317 | } | |||
318 | ||||
319 | if (!*nslices) | |||
320 | { | |||
321 | *nslices = (int)(box[axis][axis] * 10); /* default value */ | |||
322 | fprintf(stderrstderr, "\nDividing the box in %d slices\n", *nslices); | |||
323 | } | |||
324 | ||||
325 | snew(*slDensity, nr_grps)(*slDensity) = save_calloc("*slDensity", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 325, (nr_grps), sizeof(*(*slDensity))); | |||
326 | for (i = 0; i < nr_grps; i++) | |||
327 | { | |||
328 | snew((*slDensity)[i], *nslices)((*slDensity)[i]) = save_calloc("(*slDensity)[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 328, (*nslices), sizeof(*((*slDensity)[i]))); | |||
329 | } | |||
330 | ||||
331 | gpbc = gmx_rmpbc_init(&top->idef, ePBC, top->atoms.nr); | |||
332 | /*********** Start processing trajectory ***********/ | |||
333 | do | |||
334 | { | |||
335 | gmx_rmpbc(gpbc, natoms, box, x0); | |||
336 | ||||
337 | if (bCenter) | |||
338 | { | |||
339 | center_coords(&top->atoms, box, x0, axis); | |||
340 | } | |||
341 | ||||
342 | *slWidth = box[axis][axis]/(*nslices); | |||
343 | invvol = *nslices/(box[XX0][XX0]*box[YY1][YY1]*box[ZZ2][ZZ2]); | |||
344 | teller++; | |||
345 | ||||
346 | for (n = 0; n < nr_grps; n++) | |||
347 | { | |||
348 | for (i = 0; i < gnx[n]; i++) /* loop over all atoms in index file */ | |||
349 | { | |||
350 | z = x0[index[n][i]][axis]; | |||
351 | while (z < 0) | |||
352 | { | |||
353 | z += box[axis][axis]; | |||
354 | } | |||
355 | while (z > box[axis][axis]) | |||
356 | { | |||
357 | z -= box[axis][axis]; | |||
358 | } | |||
359 | ||||
360 | /* determine which slice atom is in */ | |||
361 | slice = (int)(z / (*slWidth)); | |||
362 | (*slDensity)[n][slice] += top->atoms.atom[index[n][i]].m*invvol; | |||
363 | } | |||
364 | } | |||
365 | nr_frames++; | |||
366 | } | |||
367 | while (read_next_x(oenv, status, &t, x0, box)); | |||
368 | gmx_rmpbc_done(gpbc); | |||
369 | ||||
370 | /*********** done with status file **********/ | |||
371 | close_trj(status); | |||
372 | ||||
373 | /* slDensity now contains the total mass per slice, summed over all | |||
374 | frames. Now divide by nr_frames and volume of slice | |||
375 | */ | |||
376 | ||||
377 | fprintf(stderrstderr, "\nRead %d frames from trajectory. Calculating density\n", | |||
378 | nr_frames); | |||
379 | ||||
380 | for (n = 0; n < nr_grps; n++) | |||
381 | { | |||
382 | for (i = 0; i < *nslices; i++) | |||
383 | { | |||
384 | (*slDensity)[n][i] /= nr_frames; | |||
385 | } | |||
386 | } | |||
387 | ||||
388 | sfree(x0)save_free("x0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 388, (x0)); /* free memory used by coordinate array */ | |||
389 | } | |||
390 | ||||
391 | void plot_density(double *slDensity[], const char *afile, int nslices, | |||
392 | int nr_grps, char *grpname[], real slWidth, | |||
393 | const char **dens_opt, | |||
394 | gmx_bool bSymmetrize, const output_env_t oenv) | |||
395 | { | |||
396 | FILE *den; | |||
397 | const char *ylabel = NULL((void*)0); | |||
398 | int slice, n; | |||
399 | real ddd; | |||
400 | ||||
401 | switch (dens_opt[0][0]) | |||
402 | { | |||
403 | case 'm': ylabel = "Density (kg m\\S-3\\N)"; break; | |||
404 | case 'n': ylabel = "Number density (nm\\S-3\\N)"; break; | |||
405 | case 'c': ylabel = "Charge density (e nm\\S-3\\N)"; break; | |||
406 | case 'e': ylabel = "Electron density (e nm\\S-3\\N)"; break; | |||
407 | } | |||
408 | ||||
409 | den = xvgropen(afile, "Partial densities", "Box (nm)", ylabel, oenv); | |||
410 | ||||
411 | xvgr_legend(den, nr_grps, (const char**)grpname, oenv); | |||
412 | ||||
413 | for (slice = 0; (slice < nslices); slice++) | |||
414 | { | |||
415 | fprintf(den, "%12g ", slice * slWidth); | |||
416 | for (n = 0; (n < nr_grps); n++) | |||
417 | { | |||
418 | if (bSymmetrize) | |||
419 | { | |||
420 | ddd = (slDensity[n][slice]+slDensity[n][nslices-slice-1])*0.5; | |||
421 | } | |||
422 | else | |||
423 | { | |||
424 | ddd = slDensity[n][slice]; | |||
425 | } | |||
426 | if (dens_opt[0][0] == 'm') | |||
427 | { | |||
428 | fprintf(den, " %12g", ddd*AMU(1.6605402e-27)/(NANO(1e-9)*NANO(1e-9)*NANO(1e-9))); | |||
429 | } | |||
430 | else | |||
431 | { | |||
432 | fprintf(den, " %12g", ddd); | |||
433 | } | |||
434 | } | |||
435 | fprintf(den, "\n"); | |||
436 | } | |||
437 | ||||
438 | gmx_ffclose(den); | |||
439 | } | |||
440 | ||||
441 | int gmx_density(int argc, char *argv[]) | |||
442 | { | |||
443 | const char *desc[] = { | |||
444 | "[THISMODULE] computes partial densities across the box, using an index file.[PAR]", | |||
445 | "For the total density of NPT simulations, use [gmx-energy] instead.", | |||
446 | "[PAR]", | |||
447 | "Densities are in kg/m^3, and number densities or electron densities can also be", | |||
448 | "calculated. For electron densities, a file describing the number of", | |||
449 | "electrons for each type of atom should be provided using [TT]-ei[tt].", | |||
450 | "It should look like:[BR]", | |||
451 | " [TT]2[tt][BR]", | |||
452 | " [TT]atomname = nrelectrons[tt][BR]", | |||
453 | " [TT]atomname = nrelectrons[tt][BR]", | |||
454 | "The first line contains the number of lines to read from the file.", | |||
455 | "There should be one line for each unique atom name in your system.", | |||
456 | "The number of electrons for each atom is modified by its atomic", | |||
457 | "partial charge." | |||
458 | }; | |||
459 | ||||
460 | output_env_t oenv; | |||
461 | static const char *dens_opt[] = | |||
462 | { NULL((void*)0), "mass", "number", "charge", "electron", NULL((void*)0) }; | |||
463 | static int axis = 2; /* normal to memb. default z */ | |||
464 | static const char *axtitle = "Z"; | |||
465 | static int nslices = 50; /* nr of slices defined */ | |||
466 | static int ngrps = 1; /* nr. of groups */ | |||
467 | static gmx_bool bSymmetrize = FALSE0; | |||
468 | static gmx_bool bCenter = FALSE0; | |||
469 | t_pargs pa[] = { | |||
470 | { "-d", FALSE0, etSTR, {&axtitle}, | |||
471 | "Take the normal on the membrane in direction X, Y or Z." }, | |||
472 | { "-sl", FALSE0, etINT, {&nslices}, | |||
473 | "Divide the box in this number of slices." }, | |||
474 | { "-dens", FALSE0, etENUM, {dens_opt}, | |||
475 | "Density"}, | |||
476 | { "-ng", FALSE0, etINT, {&ngrps}, | |||
477 | "Number of groups of which to compute densities." }, | |||
478 | { "-symm", FALSE0, etBOOL, {&bSymmetrize}, | |||
479 | "Symmetrize the density along the axis, with respect to the center. Useful for bilayers." }, | |||
480 | { "-center", FALSE0, etBOOL, {&bCenter}, | |||
481 | "Shift the center of mass along the axis to zero. This means if your axis is Z and your box is bX, bY, bZ, the center of mass will be at bX/2, bY/2, 0."} | |||
482 | }; | |||
483 | ||||
484 | const char *bugs[] = { | |||
485 | "When calculating electron densities, atomnames are used instead of types. This is bad.", | |||
486 | }; | |||
487 | ||||
488 | double **density; /* density per slice */ | |||
489 | real slWidth; /* width of one slice */ | |||
490 | char **grpname; /* groupnames */ | |||
491 | int nr_electrons; /* nr. electrons */ | |||
492 | int *ngx; /* sizes of groups */ | |||
493 | t_electron *el_tab; /* tabel with nr. of electrons*/ | |||
494 | t_topology *top; /* topology */ | |||
495 | int ePBC; | |||
496 | atom_id **index; /* indices for all groups */ | |||
497 | int i; | |||
498 | ||||
499 | t_filenm fnm[] = { /* files for g_density */ | |||
500 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, | |||
501 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
502 | { efTPX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, | |||
503 | { efDAT, "-ei", "electrons", ffOPTRD(1<<1 | 1<<3) }, /* file with nr. of electrons */ | |||
504 | { efXVG, "-o", "density", ffWRITE1<<2 }, | |||
505 | }; | |||
506 | ||||
507 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) | |||
508 | ||||
509 | 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), | |||
| ||||
510 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, asize(bugs)((int)(sizeof(bugs)/sizeof((bugs)[0]))), bugs, | |||
511 | &oenv)) | |||
512 | { | |||
513 | return 0; | |||
514 | } | |||
515 | ||||
516 | if (bSymmetrize && !bCenter) | |||
517 | { | |||
518 | fprintf(stderrstderr, "Can not symmetrize without centering. Turning on -center\n"); | |||
519 | bCenter = TRUE1; | |||
520 | } | |||
521 | /* Calculate axis */ | |||
522 | axis = toupper(axtitle[0])(__extension__ ({ int __res; if (sizeof (axtitle[0]) > 1) { if (__builtin_constant_p (axtitle[0])) { int __c = (axtitle[ 0]); __res = __c < -128 || __c > 255 ? __c : (*__ctype_toupper_loc ())[__c]; } else __res = toupper (axtitle[0]); } else __res = (*__ctype_toupper_loc ())[(int) (axtitle[0])]; __res; })) - 'X'; | |||
523 | ||||
524 | top = read_top(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &ePBC); /* read topology file */ | |||
525 | if (dens_opt[0][0] == 'n') | |||
| ||||
526 | { | |||
527 | for (i = 0; (i < top->atoms.nr); i++) | |||
528 | { | |||
529 | top->atoms.atom[i].m = 1; | |||
530 | } | |||
531 | } | |||
532 | else if (dens_opt[0][0] == 'c') | |||
533 | { | |||
534 | for (i = 0; (i < top->atoms.nr); i++) | |||
535 | { | |||
536 | top->atoms.atom[i].m = top->atoms.atom[i].q; | |||
537 | } | |||
538 | } | |||
539 | ||||
540 | snew(grpname, ngrps)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 540, (ngrps), sizeof(*(grpname))); | |||
541 | snew(index, ngrps)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 541, (ngrps), sizeof(*(index))); | |||
542 | snew(ngx, ngrps)(ngx) = save_calloc("ngx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_density.c" , 542, (ngrps), sizeof(*(ngx))); | |||
543 | ||||
544 | get_index(&top->atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), ngrps, ngx, index, grpname); | |||
545 | ||||
546 | if (dens_opt[0][0] == 'e') | |||
547 | { | |||
548 | nr_electrons = get_electrons(&el_tab, ftp2fn(efDAT, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); | |||
549 | fprintf(stderrstderr, "Read %d atomtypes from datafile\n", nr_electrons); | |||
550 | ||||
551 | calc_electron_density(ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), index, ngx, &density, | |||
552 | &nslices, top, ePBC, axis, ngrps, &slWidth, el_tab, | |||
553 | nr_electrons, bCenter, oenv); | |||
554 | } | |||
555 | else | |||
556 | { | |||
557 | calc_density(ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), index, ngx, &density, &nslices, top, | |||
558 | ePBC, axis, ngrps, &slWidth, bCenter, oenv); | |||
559 | } | |||
560 | ||||
561 | plot_density(density, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
562 | nslices, ngrps, grpname, slWidth, dens_opt, | |||
563 | bSymmetrize, oenv); | |||
564 | ||||
565 | do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy"); /* view xvgr file */ | |||
566 | return 0; | |||
567 | } |