Bug Summary

File:gromacs/gmxana/gmx_mindist.c
Location:line 580, column 5
Description:Access to field 'nres' results in a dereference of a null pointer (loaded from variable 'atoms')

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 <stdlib.h>
43#include <string.h>
44
45#include "typedefs.h"
46#include "gromacs/utility/smalloc.h"
47#include "macros.h"
48#include "gromacs/math/vec.h"
49#include "gromacs/fileio/xvgr.h"
50#include "viewit.h"
51#include "pbc.h"
52#include "gromacs/utility/futil.h"
53#include "gromacs/commandline/pargs.h"
54#include "index.h"
55#include "gromacs/fileio/tpxio.h"
56#include "gromacs/fileio/trxio.h"
57#include "rmpbc.h"
58#include "gmx_ana.h"
59
60
61static void periodic_dist(matrix box, rvec x[], int n, atom_id index[],
62 real *rmin, real *rmax, int *min_ind)
63{
64#define NSHIFT26 26
65 int sx, sy, sz, i, j, s;
66 real sqr_box, r2min, r2max, r2;
67 rvec shift[NSHIFT26], d0, d;
68
69 sqr_box = sqr(min(norm(box[XX]), min(norm(box[YY]), norm(box[ZZ])))(((norm(box[0])) < ((((norm(box[1])) < (norm(box[2]))) ?
(norm(box[1])) : (norm(box[2])) ))) ? (norm(box[0])) : ((((norm
(box[1])) < (norm(box[2]))) ? (norm(box[1])) : (norm(box[2
])) )) ));
70
71 s = 0;
72 for (sz = -1; sz <= 1; sz++)
73 {
74 for (sy = -1; sy <= 1; sy++)
75 {
76 for (sx = -1; sx <= 1; sx++)
77 {
78 if (sx != 0 || sy != 0 || sz != 0)
79 {
80 for (i = 0; i < DIM3; i++)
81 {
82 shift[s][i] = sx*box[XX0][i]+sy*box[YY1][i]+sz*box[ZZ2][i];
83 }
84 s++;
85 }
86 }
87 }
88 }
89
90 r2min = sqr_box;
91 r2max = 0;
92
93 for (i = 0; i < n; i++)
94 {
95 for (j = i+1; j < n; j++)
96 {
97 rvec_sub(x[index[i]], x[index[j]], d0);
98 r2 = norm2(d0);
99 if (r2 > r2max)
100 {
101 r2max = r2;
102 }
103 for (s = 0; s < NSHIFT26; s++)
104 {
105 rvec_add(d0, shift[s], d);
106 r2 = norm2(d);
107 if (r2 < r2min)
108 {
109 r2min = r2;
110 min_ind[0] = i;
111 min_ind[1] = j;
112 }
113 }
114 }
115 }
116
117 *rmin = sqrt(r2min);
118 *rmax = sqrt(r2max);
119}
120
121static void periodic_mindist_plot(const char *trxfn, const char *outfn,
122 t_topology *top, int ePBC,
123 int n, atom_id index[], gmx_bool bSplit,
124 const output_env_t oenv)
125{
126 FILE *out;
127 const char *leg[5] = { "min per.", "max int.", "box1", "box2", "box3" };
128 t_trxstatus *status;
129 real t;
130 rvec *x;
131 matrix box;
132 int natoms, ind_min[2] = {0, 0}, ind_mini = 0, ind_minj = 0;
133 real r, rmin, rmax, rmint, tmint;
134 gmx_bool bFirst;
135 gmx_rmpbc_t gpbc = NULL((void*)0);
136
137 natoms = read_first_x(oenv, &status, trxfn, &t, &x, box);
138
139 check_index(NULL((void*)0), n, index, NULL((void*)0), natoms);
140
141 out = xvgropen(outfn, "Minimum distance to periodic image",
142 output_env_get_time_label(oenv), "Distance (nm)", oenv);
143 if (output_env_get_print_xvgr_codes(oenv))
144 {
145 fprintf(out, "@ subtitle \"and maximum internal distance\"\n");
146 }
147 xvgr_legend(out, 5, leg, oenv);
148
149 rmint = box[XX0][XX0];
150 tmint = 0;
151
152 if (NULL((void*)0) != top)
153 {
154 gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
155 }
156
157 bFirst = TRUE1;
158 do
159 {
160 if (NULL((void*)0) != top)
161 {
162 gmx_rmpbc(gpbc, natoms, box, x);
163 }
164
165 periodic_dist(box, x, n, index, &rmin, &rmax, ind_min);
166 if (rmin < rmint)
167 {
168 rmint = rmin;
169 tmint = t;
170 ind_mini = ind_min[0];
171 ind_minj = ind_min[1];
172 }
173 if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
174 {
175 fprintf(out, "&\n");
176 }
177 fprintf(out, "\t%g\t%6.3f %6.3f %6.3f %6.3f %6.3f\n",
178 output_env_conv_time(oenv, t), rmin, rmax, norm(box[0]), norm(box[1]), norm(box[2]));
179 bFirst = FALSE0;
180 }
181 while (read_next_x(oenv, status, &t, x, box));
182
183 if (NULL((void*)0) != top)
184 {
185 gmx_rmpbc_done(gpbc);
186 }
187
188 gmx_ffclose(out);
189
190 fprintf(stdoutstdout,
191 "\nThe shortest periodic distance is %g (nm) at time %g (%s),\n"
192 "between atoms %d and %d\n",
193 rmint, output_env_conv_time(oenv, tmint), output_env_get_time_unit(oenv),
194 index[ind_mini]+1, index[ind_minj]+1);
195}
196
197static void calc_dist(real rcut, gmx_bool bPBC, int ePBC, matrix box, rvec x[],
198 int nx1, int nx2, atom_id index1[], atom_id index2[],
199 gmx_bool bGroup,
200 real *rmin, real *rmax, int *nmin, int *nmax,
201 int *ixmin, int *jxmin, int *ixmax, int *jxmax)
202{
203 int i, j, i0 = 0, j1;
204 int ix, jx;
205 atom_id *index3;
206 rvec dx;
207 real r2, rmin2, rmax2, rcut2;
208 t_pbc pbc;
209 int nmin_j, nmax_j;
210
211 *ixmin = -1;
212 *jxmin = -1;
213 *ixmax = -1;
214 *jxmax = -1;
215 *nmin = 0;
216 *nmax = 0;
217
218 rcut2 = sqr(rcut);
219
220 /* Must init pbc every step because of pressure coupling */
221 if (bPBC)
222 {
223 set_pbc(&pbc, ePBC, box);
224 }
225 if (index2)
226 {
227 i0 = 0;
228 j1 = nx2;
229 index3 = index2;
230 }
231 else
232 {
233 j1 = nx1;
234 index3 = index1;
235 }
236
237 rmin2 = 1e12;
238 rmax2 = -1e12;
239
240 for (j = 0; (j < j1); j++)
241 {
242 jx = index3[j];
243 if (index2 == NULL((void*)0))
244 {
245 i0 = j + 1;
246 }
247 nmin_j = 0;
248 nmax_j = 0;
249 for (i = i0; (i < nx1); i++)
250 {
251 ix = index1[i];
252 if (ix != jx)
253 {
254 if (bPBC)
255 {
256 pbc_dx(&pbc, x[ix], x[jx], dx);
257 }
258 else
259 {
260 rvec_sub(x[ix], x[jx], dx);
261 }
262 r2 = iprod(dx, dx);
263 if (r2 < rmin2)
264 {
265 rmin2 = r2;
266 *ixmin = ix;
267 *jxmin = jx;
268 }
269 if (r2 > rmax2)
270 {
271 rmax2 = r2;
272 *ixmax = ix;
273 *jxmax = jx;
274 }
275 if (r2 <= rcut2)
276 {
277 nmin_j++;
278 }
279 else if (r2 > rcut2)
280 {
281 nmax_j++;
282 }
283 }
284 }
285 if (bGroup)
286 {
287 if (nmin_j > 0)
288 {
289 (*nmin)++;
290 }
291 if (nmax_j > 0)
292 {
293 (*nmax)++;
294 }
295 }
296 else
297 {
298 *nmin += nmin_j;
299 *nmax += nmax_j;
300 }
301 }
302 *rmin = sqrt(rmin2);
303 *rmax = sqrt(rmax2);
304}
305
306void dist_plot(const char *fn, const char *afile, const char *dfile,
307 const char *nfile, const char *rfile, const char *xfile,
308 real rcut, gmx_bool bMat, t_atoms *atoms,
309 int ng, atom_id *index[], int gnx[], char *grpn[], gmx_bool bSplit,
310 gmx_bool bMin, int nres, atom_id *residue, gmx_bool bPBC, int ePBC,
311 gmx_bool bGroup, gmx_bool bEachResEachTime, gmx_bool bPrintResName,
312 const output_env_t oenv)
313{
314 FILE *atm, *dist, *num;
315 t_trxstatus *trxout;
316 char buf[256];
317 char **leg;
318 real t, dmin, dmax, **mindres = NULL((void*)0), **maxdres = NULL((void*)0);
319 int nmin, nmax;
320 t_trxstatus *status;
321 int i = -1, j, k, natoms;
322 int min1, min2, max1, max2, min1r, min2r, max1r, max2r;
323 atom_id oindex[2];
324 rvec *x0;
325 matrix box;
326 t_trxframe frout;
327 gmx_bool bFirst;
328 FILE *respertime = NULL((void*)0);
329
330 if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
331 {
332 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 332, "Could not read coordinates from statusfile\n");
333 }
334
335 sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
336 dist = xvgropen(dfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
337 sprintf(buf, "Number of Contacts %s %g nm", bMin ? "<" : ">", rcut);
338 num = nfile ? xvgropen(nfile, buf, output_env_get_time_label(oenv), "Number", oenv) : NULL((void*)0);
339 atm = afile ? gmx_ffopen(afile, "w") : NULL((void*)0);
340 trxout = xfile ? open_trx(xfile, "w") : NULL((void*)0);
341
342 if (bMat)
343 {
344 if (ng == 1)
345 {
346 snew(leg, 1)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 346, (1), sizeof(*(leg)));
347 sprintf(buf, "Internal in %s", grpn[0]);
348 leg[0] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t
)(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1
) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, buf, __len); __retval
; })) : __strdup (buf)));
349 xvgr_legend(dist, 0, (const char**)leg, oenv);
350 if (num)
351 {
352 xvgr_legend(num, 0, (const char**)leg, oenv);
353 }
354 }
355 else
356 {
357 snew(leg, (ng*(ng-1))/2)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 357, ((ng*(ng-1))/2), sizeof(*(leg)));
358 for (i = j = 0; (i < ng-1); i++)
359 {
360 for (k = i+1; (k < ng); k++, j++)
361 {
362 sprintf(buf, "%s-%s", grpn[i], grpn[k]);
363 leg[j] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t
)(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1
) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, buf, __len); __retval
; })) : __strdup (buf)));
364 }
365 }
366 xvgr_legend(dist, j, (const char**)leg, oenv);
367 if (num)
368 {
369 xvgr_legend(num, j, (const char**)leg, oenv);
370 }
371 }
372 }
373 else
374 {
375 snew(leg, ng-1)(leg) = save_calloc("leg", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 375, (ng-1), sizeof(*(leg)));
376 for (i = 0; (i < ng-1); i++)
377 {
378 sprintf(buf, "%s-%s", grpn[0], grpn[i+1]);
379 leg[i] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t
)(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1
) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, buf, __len); __retval
; })) : __strdup (buf)));
380 }
381 xvgr_legend(dist, ng-1, (const char**)leg, oenv);
382 if (num)
383 {
384 xvgr_legend(num, ng-1, (const char**)leg, oenv);
385 }
386 }
387
388 if (bEachResEachTime)
389 {
390 sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
391 respertime = xvgropen(rfile, buf, output_env_get_time_label(oenv), "Distance (nm)", oenv);
392 xvgr_legend(respertime, ng-1, (const char**)leg, oenv);
393 if (bPrintResName)
394 {
395 fprintf(respertime, "# ");
396 }
397 for (j = 0; j < nres; j++)
398 {
399 fprintf(respertime, "%s%d ", *(atoms->resinfo[atoms->atom[index[0][residue[j]]].resind].name), atoms->atom[index[0][residue[j]]].resind);
400 }
401 fprintf(respertime, "\n");
402
403 }
404
405 j = 0;
406 if (nres)
407 {
408 snew(mindres, ng-1)(mindres) = save_calloc("mindres", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 408, (ng-1), sizeof(*(mindres)));
409 snew(maxdres, ng-1)(maxdres) = save_calloc("maxdres", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 409, (ng-1), sizeof(*(maxdres)));
410 for (i = 1; i < ng; i++)
411 {
412 snew(mindres[i-1], nres)(mindres[i-1]) = save_calloc("mindres[i-1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 412, (nres), sizeof(*(mindres[i-1])));
413 snew(maxdres[i-1], nres)(maxdres[i-1]) = save_calloc("maxdres[i-1]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 413, (nres), sizeof(*(maxdres[i-1])));
414 for (j = 0; j < nres; j++)
415 {
416 mindres[i-1][j] = 1e6;
417 }
418 /* maxdres[*][*] is already 0 */
419 }
420 }
421 bFirst = TRUE1;
422 do
423 {
424 if (bSplit && !bFirst && abs(t/output_env_get_time_factor(oenv)) < 1e-5)
425 {
426 fprintf(dist, "&\n");
427 if (num)
428 {
429 fprintf(num, "&\n");
430 }
431 if (atm)
432 {
433 fprintf(atm, "&\n");
434 }
435 }
436 fprintf(dist, "%12e", output_env_conv_time(oenv, t));
437 if (num)
438 {
439 fprintf(num, "%12e", output_env_conv_time(oenv, t));
440 }
441
442 if (bMat)
443 {
444 if (ng == 1)
445 {
446 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[0], index[0], index[0], bGroup,
447 &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
448 fprintf(dist, " %12e", bMin ? dmin : dmax);
449 if (num)
450 {
451 fprintf(num, " %8d", bMin ? nmin : nmax);
452 }
453 }
454 else
455 {
456 for (i = 0; (i < ng-1); i++)
457 {
458 for (k = i+1; (k < ng); k++)
459 {
460 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[i], gnx[k], index[i], index[k],
461 bGroup, &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
462 fprintf(dist, " %12e", bMin ? dmin : dmax);
463 if (num)
464 {
465 fprintf(num, " %8d", bMin ? nmin : nmax);
466 }
467 }
468 }
469 }
470 }
471 else
472 {
473 for (i = 1; (i < ng); i++)
474 {
475 calc_dist(rcut, bPBC, ePBC, box, x0, gnx[0], gnx[i], index[0], index[i], bGroup,
476 &dmin, &dmax, &nmin, &nmax, &min1, &min2, &max1, &max2);
477 fprintf(dist, " %12e", bMin ? dmin : dmax);
478 if (num)
479 {
480 fprintf(num, " %8d", bMin ? nmin : nmax);
481 }
482 if (nres)
483 {
484 for (j = 0; j < nres; j++)
485 {
486 calc_dist(rcut, bPBC, ePBC, box, x0, residue[j+1]-residue[j], gnx[i],
487 &(index[0][residue[j]]), index[i], bGroup,
488 &dmin, &dmax, &nmin, &nmax, &min1r, &min2r, &max1r, &max2r);
489 mindres[i-1][j] = min(mindres[i-1][j], dmin)(((mindres[i-1][j]) < (dmin)) ? (mindres[i-1][j]) : (dmin)
);
490 maxdres[i-1][j] = max(maxdres[i-1][j], dmax)(((maxdres[i-1][j]) > (dmax)) ? (maxdres[i-1][j]) : (dmax)
);
491 }
492 }
493 }
494 }
495 fprintf(dist, "\n");
496 if (num)
497 {
498 fprintf(num, "\n");
499 }
500 if ( (bMin ? min1 : max1) != -1)
501 {
502 if (atm)
503 {
504 fprintf(atm, "%12e %12d %12d\n",
505 output_env_conv_time(oenv, t), 1+(bMin ? min1 : max1),
506 1+(bMin ? min2 : max2));
507 }
508 }
509
510 if (trxout)
511 {
512 oindex[0] = bMin ? min1 : max1;
513 oindex[1] = bMin ? min2 : max2;
514 write_trx(trxout, 2, oindex, atoms, i, t, box, x0, NULL((void*)0), NULL((void*)0));
515 }
516 bFirst = FALSE0;
517 /*dmin should be minimum distance for residue and group*/
518 if (bEachResEachTime)
519 {
520 fprintf(respertime, "%12e", t);
521 for (i = 1; i < ng; i++)
522 {
523 for (j = 0; j < nres; j++)
524 {
525 fprintf(respertime, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
526 /*reset distances for next time point*/
527 mindres[i-1][j] = 1e6;
528 maxdres[i-1][j] = 0;
529 }
530 }
531 fprintf(respertime, "\n");
532 }
533 }
534 while (read_next_x(oenv, status, &t, x0, box));
535
536 close_trj(status);
537 gmx_ffclose(dist);
538 if (num)
539 {
540 gmx_ffclose(num);
541 }
542 if (atm)
543 {
544 gmx_ffclose(atm);
545 }
546 if (trxout)
547 {
548 close_trx(trxout);
549 }
550
551 if (nres && !bEachResEachTime)
552 {
553 FILE *res;
554
555 sprintf(buf, "%simum Distance", bMin ? "Min" : "Max");
556 res = xvgropen(rfile, buf, "Residue (#)", "Distance (nm)", oenv);
557 xvgr_legend(res, ng-1, (const char**)leg, oenv);
558 for (j = 0; j < nres; j++)
559 {
560 fprintf(res, "%4d", j+1);
561 for (i = 1; i < ng; i++)
562 {
563 fprintf(res, " %7g", bMin ? mindres[i-1][j] : maxdres[i-1][j]);
564 }
565 fprintf(res, "\n");
566 }
567 }
568
569 sfree(x0)save_free("x0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 569, (x0));
570}
571
572int find_residues(t_atoms *atoms, int n, atom_id index[], atom_id **resindex)
573{
574 int i;
575 int nres = 0, resnr, presnr;
576 int *residx;
577
578 /* build index of first atom numbers for each residue */
579 presnr = NOTSET-12345;
580 snew(residx, atoms->nres+1)(residx) = save_calloc("residx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 580, (atoms->nres+1), sizeof(*(residx)));
13
←
Within the expansion of the macro 'snew':
a
Access to field 'nres' results in a dereference of a null pointer (loaded from variable 'atoms')
581 for (i = 0; i < n; i++)
582 {
583 resnr = atoms->atom[index[i]].resind;
584 if (resnr != presnr)
585 {
586 residx[nres] = i;
587 nres++;
588 presnr = resnr;
589 }
590 }
591 if (debug)
592 {
593 printf("Found %d residues out of %d (%d/%d atoms)\n",
594 nres, atoms->nres, atoms->nr, n);
595 }
596 srenew(residx, nres+1)(residx) = save_realloc("residx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 596, (residx), (nres+1), sizeof(*(residx)));
597 /* mark end of last residue */
598 residx[nres] = n;
599 *resindex = residx;
600 return nres;
601}
602
603void dump_res(FILE *out, int nres, atom_id *resindex, atom_id index[])
604{
605 int i, j;
606
607 for (i = 0; i < nres-1; i++)
608 {
609 fprintf(out, "Res %d (%d):", i, resindex[i+1]-resindex[i]);
610 for (j = resindex[i]; j < resindex[i+1]; j++)
611 {
612 fprintf(out, " %d(%d)", j, index[j]);
613 }
614 fprintf(out, "\n");
615 }
616}
617
618int gmx_mindist(int argc, char *argv[])
619{
620 const char *desc[] = {
621 "[THISMODULE] computes the distance between one group and a number of",
622 "other groups. Both the minimum distance",
623 "(between any pair of atoms from the respective groups)",
624 "and the number of contacts within a given",
625 "distance are written to two separate output files.",
626 "With the [TT]-group[tt] option a contact of an atom in another group",
627 "with multiple atoms in the first group is counted as one contact",
628 "instead of as multiple contacts.",
629 "With [TT]-or[tt], minimum distances to each residue in the first",
630 "group are determined and plotted as a function of residue number.[PAR]",
631 "With option [TT]-pi[tt] the minimum distance of a group to its",
632 "periodic image is plotted. This is useful for checking if a protein",
633 "has seen its periodic image during a simulation. Only one shift in",
634 "each direction is considered, giving a total of 26 shifts.",
635 "It also plots the maximum distance within the group and the lengths",
636 "of the three box vectors.[PAR]",
637 "Also [gmx-distance] calculates distances."
638 };
639 const char *bugs[] = {
640 "The [TT]-pi[tt] option is very slow."
641 };
642
643 static gmx_bool bMat = FALSE0, bPI = FALSE0, bSplit = FALSE0, bMax = FALSE0, bPBC = TRUE1;
644 static gmx_bool bGroup = FALSE0;
645 static real rcutoff = 0.6;
646 static int ng = 1;
647 static gmx_bool bEachResEachTime = FALSE0, bPrintResName = FALSE0;
648 t_pargs pa[] = {
649 { "-matrix", FALSE0, etBOOL, {&bMat},
650 "Calculate half a matrix of group-group distances" },
651 { "-max", FALSE0, etBOOL, {&bMax},
652 "Calculate *maximum* distance instead of minimum" },
653 { "-d", FALSE0, etREAL, {&rcutoff},
654 "Distance for contacts" },
655 { "-group", FALSE0, etBOOL, {&bGroup},
656 "Count contacts with multiple atoms in the first group as one" },
657 { "-pi", FALSE0, etBOOL, {&bPI},
658 "Calculate minimum distance with periodic images" },
659 { "-split", FALSE0, etBOOL, {&bSplit},
660 "Split graph where time is zero" },
661 { "-ng", FALSE0, etINT, {&ng},
662 "Number of secondary groups to compute distance to a central group" },
663 { "-pbc", FALSE0, etBOOL, {&bPBC},
664 "Take periodic boundary conditions into account" },
665 { "-respertime", FALSE0, etBOOL, {&bEachResEachTime},
666 "When writing per-residue distances, write distance for each time point" },
667 { "-printresname", FALSE0, etBOOL, {&bPrintResName},
668 "Write residue names" }
669 };
670 output_env_t oenv;
671 t_topology *top = NULL((void*)0);
672 int ePBC = -1;
673 char title[256];
674 real t;
675 rvec *x;
676 matrix box;
677 gmx_bool bTop = FALSE0;
678
679 FILE *atm;
680 int i, j, nres = 0;
681 const char *trxfnm, *tpsfnm, *ndxfnm, *distfnm, *numfnm, *atmfnm, *oxfnm, *resfnm;
682 char **grpname;
683 int *gnx;
684 atom_id **index, *residues = NULL((void*)0);
685 t_filenm fnm[] = {
686 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
687 { efTPS, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
688 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
689 { efXVG, "-od", "mindist", ffWRITE1<<2 },
690 { efXVG, "-on", "numcont", ffOPTWR(1<<2| 1<<3) },
691 { efOUT, "-o", "atm-pair", ffOPTWR(1<<2| 1<<3) },
692 { efTRO, "-ox", "mindist", ffOPTWR(1<<2| 1<<3) },
693 { efXVG, "-or", "mindistres", ffOPTWR(1<<2| 1<<3) }
694 };
695#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
696
697 if (!parse_common_args(&argc, argv,
1
Taking false branch
→
698 PCA_CAN_VIEW(1<<5) | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13),
699 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, 0, NULL((void*)0), &oenv))
700 {
701 return 0;
702 }
703
704 trxfnm = ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
705 ndxfnm = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
706 distfnm = opt2fn("-od", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
707 numfnm = opt2fn_null("-on", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
708 atmfnm = ftp2fn_null(efOUT, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
709 oxfnm = opt2fn_null("-ox", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
710 resfnm = opt2fn_null("-or", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
711 if (bPI || resfnm != NULL((void*)0))
2
←
Assuming 'resfnm' is not equal to null
→
3
←
Taking true branch
→
712 {
713 /* We need a tps file */
714 tpsfnm = ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
715 }
716 else
717 {
718 tpsfnm = ftp2fn_null(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
719 }
720
721 if (!tpsfnm && !ndxfnm)
4
←
Assuming 'tpsfnm' is non-null
→
722 {
723 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 723, "You have to specify either the index file or a tpr file");
724 }
725
726 if (bPI)
5
←
Taking false branch
→
727 {
728 ng = 1;
729 fprintf(stderrstderr, "Choose a group for distance calculation\n");
730 }
731 else if (!bMat)
6
←
Taking true branch
→
732 {
733 ng++;
734 }
735
736 snew(gnx, ng)(gnx) = save_calloc("gnx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 736, (ng), sizeof(*(gnx)));
737 snew(index, ng)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 737, (ng), sizeof(*(index)));
738 snew(grpname, ng)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 738, (ng), sizeof(*(grpname)));
739
740 if (tpsfnm || resfnm || !ndxfnm)
741 {
742 snew(top, 1)(top) = save_calloc("top", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 742, (1), sizeof(*(top)));
743 bTop = read_tps_conf(tpsfnm, title, top, &ePBC, &x, NULL((void*)0), box, FALSE0);
744 if (bPI && !bTop)
745 {
746 printf("\nWARNING: Without a run input file a trajectory with broken molecules will not give the correct periodic image distance\n\n");
747 }
748 }
749 get_index(top ? &(top->atoms) : NULL((void*)0), ndxfnm, ng, gnx, index, grpname);
7
←
Assuming 'top' is null
→
8
←
'?' condition is false
→
750
751 if (bMat && (ng == 1))
752 {
753 ng = gnx[0];
754 printf("Special case: making distance matrix between all atoms in group %s\n",
755 grpname[0]);
756 srenew(gnx, ng)(gnx) = save_realloc("gnx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 756, (gnx), (ng), sizeof(*(gnx)));
757 srenew(index, ng)(index) = save_realloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 757, (index), (ng), sizeof(*(index)));
758 srenew(grpname, ng)(grpname) = save_realloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 758, (grpname), (ng), sizeof(*(grpname)));
759 for (i = 1; (i < ng); i++)
760 {
761 gnx[i] = 1;
762 grpname[i] = grpname[0];
763 snew(index[i], 1)(index[i]) = save_calloc("index[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_mindist.c"
, 763, (1), sizeof(*(index[i])));
764 index[i][0] = index[0][i];
765 }
766 gnx[0] = 1;
767 }
768
769 if (resfnm)
9
←
Taking true branch
→
770 {
771 nres = find_residues(top ? &(top->atoms) : NULL((void*)0),
10
←
'?' condition is false
→
11
←
Passing null pointer value via 1st parameter 'atoms'
→
12
←
Calling 'find_residues'
→
772 gnx[0], index[0], &residues);
773 if (debug)
774 {
775 dump_res(debug, nres, residues, index[0]);
776 }
777 }
778
779 if (bPI)
780 {
781 periodic_mindist_plot(trxfnm, distfnm, top, ePBC, gnx[0], index[0], bSplit, oenv);
782 }
783 else
784 {
785 dist_plot(trxfnm, atmfnm, distfnm, numfnm, resfnm, oxfnm,
786 rcutoff, bMat, top ? &(top->atoms) : NULL((void*)0),
787 ng, index, gnx, grpname, bSplit, !bMax, nres, residues, bPBC, ePBC,
788 bGroup, bEachResEachTime, bPrintResName, oenv);
789 }
790
791 do_view(oenv, distfnm, "-nxy");
792 if (!bPI)
793 {
794 do_view(oenv, numfnm, "-nxy");
795 }
796
797 return 0;
798}