File: | gromacs/gmxana/gmx_mindist.c |
Location: | line 500, column 35 |
Description: | The left operand of '!=' is a garbage value |
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[] = { | |||
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 | } |