File: | gromacs/gmxana/gmx_mindist.c |
Location: | line 639, column 21 |
Description: | Value stored to 'bugs' during its initialization 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 <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 | |
61 | static 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 | |
121 | static 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 | |
197 | static 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 | |
306 | void 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 | |
572 | int 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))); |
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 | |
603 | void 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 | |
618 | int 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[] = { |
Value stored to 'bugs' during its initialization is never read | |
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, |
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)) |
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) |
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) |
727 | { |
728 | ng = 1; |
729 | fprintf(stderrstderr, "Choose a group for distance calculation\n"); |
730 | } |
731 | else if (!bMat) |
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); |
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) |
770 | { |
771 | nres = find_residues(top ? &(top->atoms) : NULL((void*)0), |
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 | } |