Bug Summary

File:gromacs/gmxana/gmx_confrms.c
Location:line 136, column 5
Description:Value stored to 'dx' is never read

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 <string.h>
43
44#include "gromacs/fileio/filenm.h"
45#include "gromacs/utility/smalloc.h"
46#include "macros.h"
47#include "typedefs.h"
48#include "gromacs/commandline/pargs.h"
49#include "gromacs/fileio/tpxio.h"
50#include "gromacs/math/vec.h"
51#include "index.h"
52#include "pbc.h"
53#include "gromacs/utility/fatalerror.h"
54#include "gromacs/utility/futil.h"
55#include "gromacs/fileio/confio.h"
56#include "gromacs/fileio/pdbio.h"
57#include "txtdump.h"
58#include "viewit.h"
59#include "rmpbc.h"
60#include "gmx_ana.h"
61
62#include "gromacs/math/do_fit.h"
63
64void calc_rm_cm(int isize, atom_id index[], t_atoms *atoms, rvec x[], rvec xcm)
65{
66 int i, d;
67 real tm, m;
68
69 /* calculate and remove center of mass of reference structure */
70 tm = 0;
71 clear_rvec(xcm);
72 for (i = 0; i < isize; i++)
73 {
74 m = atoms->atom[index[i]].m;
75 for (d = 0; d < DIM3; d++)
76 {
77 xcm[d] += m*x[index[i]][d];
78 }
79 tm += m;
80 }
81 svmul(1/tm, xcm, xcm);
82 for (i = 0; i < atoms->nr; i++)
83 {
84 rvec_dec(x[i], xcm);
85 }
86}
87
88int build_res_index(int isize, atom_id index[], t_atom atom[], int rindex[])
89{
90 int i, r;
91
92 r = 0;
93 rindex[r] = atom[index[0]].resind;
94 r++;
95 for (i = 1; i < isize; i++)
96 {
97 if (atom[index[i]].resind != rindex[r-1])
98 {
99 rindex[r] = atom[index[i]].resind;
100 r++;
101 }
102 }
103
104 return r;
105}
106
107int find_res_end(int i, int isize, atom_id index[], t_atoms *atoms)
108{
109 int rnr;
110
111 rnr = atoms->atom[index[i]].resind;
112 while (i < isize && atoms->atom[index[i]].resind == rnr)
113 {
114 i++;
115 }
116 return i;
117}
118
119int debug_strcmp(char s1[], char s2[])
120{
121 if (debug)
122 {
123 fprintf(debug, " %s-%s", s1, s2);
124 }
125 return strcmp(s1, s2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(s1) && __builtin_constant_p (s2) && (__s1_len
= strlen (s1), __s2_len = strlen (s2), (!((size_t)(const void
*)((s1) + 1) - (size_t)(const void *)(s1) == 1) || __s1_len >=
4) && (!((size_t)(const void *)((s2) + 1) - (size_t)
(const void *)(s2) == 1) || __s2_len >= 4)) ? __builtin_strcmp
(s1, s2) : (__builtin_constant_p (s1) && ((size_t)(const
void *)((s1) + 1) - (size_t)(const void *)(s1) == 1) &&
(__s1_len = strlen (s1), __s1_len < 4) ? (__builtin_constant_p
(s2) && ((size_t)(const void *)((s2) + 1) - (size_t)
(const void *)(s2) == 1) ? __builtin_strcmp (s1, s2) : (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (s2); int __result = (((const unsigned char *) (const
char *) (s1))[0] - __s2[0]); if (__s1_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
s1))[1] - __s2[1]); if (__s1_len > 1 && __result ==
0) { __result = (((const unsigned char *) (const char *) (s1
))[2] - __s2[2]); if (__s1_len > 2 && __result == 0
) __result = (((const unsigned char *) (const char *) (s1))[3
] - __s2[3]); } } __result; }))) : (__builtin_constant_p (s2)
&& ((size_t)(const void *)((s2) + 1) - (size_t)(const
void *)(s2) == 1) && (__s2_len = strlen (s2), __s2_len
< 4) ? (__builtin_constant_p (s1) && ((size_t)(const
void *)((s1) + 1) - (size_t)(const void *)(s1) == 1) ? __builtin_strcmp
(s1, s2) : (- (__extension__ ({ const unsigned char *__s2 = (
const unsigned char *) (const char *) (s1); int __result = ((
(const unsigned char *) (const char *) (s2))[0] - __s2[0]); if
(__s2_len > 0 && __result == 0) { __result = (((const
unsigned char *) (const char *) (s2))[1] - __s2[1]); if (__s2_len
> 1 && __result == 0) { __result = (((const unsigned
char *) (const char *) (s2))[2] - __s2[2]); if (__s2_len >
2 && __result == 0) __result = (((const unsigned char
*) (const char *) (s2))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (s1, s2)))); });
126}
127
128int find_next_match_atoms_in_res(int *i1, atom_id index1[],
129 int m1, char **atnms1[],
130 int *i2, atom_id index2[],
131 int m2, char **atnms2[])
132{
133 int dx, dy, dmax, cmp;
134 gmx_bool bFW = FALSE0;
135
136 dx = dy = 0;
Value stored to 'dx' is never read
137 cmp = NOTSET-12345;
138 dmax = max(m1-*i1, m2-*i2)(((m1-*i1) > (m2-*i2)) ? (m1-*i1) : (m2-*i2) );
139 for (dx = 0; dx < dmax && cmp != 0; dx++)
140 {
141 for (dy = dx; dy < dmax && cmp != 0; dy++)
142 {
143 if (dx || dy)
144 {
145 if (debug)
146 {
147 fprintf(debug, ".");
148 }
149 cmp = NOTSET-12345;
150 if (*i1+dx < m1 && *i2+dy < m2)
151 {
152 bFW = TRUE1;
153 cmp = debug_strcmp(*atnms1[index1[*i1+dx]], *atnms2[index2[*i2+dy]]);
154 if (debug)
155 {
156 fprintf(debug, "(%d %d)", *i1+dx, *i2+dy);
157 }
158 }
159 if (cmp != 0 && *i1+dy < m1 && *i2+dx < m2)
160 {
161 bFW = FALSE0;
162 cmp = debug_strcmp(*atnms1[index1[*i1+dy]], *atnms2[index2[*i2+dx]]);
163 if (debug)
164 {
165 fprintf(debug, "(%d %d)", *i1+dy, *i2+dx);
166 }
167 }
168 }
169 }
170 }
171 /* apparently, dx and dy are incremented one more time
172 as the loops terminate: we correct this here */
173 dx--;
174 dy--;
175 if (cmp == 0)
176 {
177 if (debug)
178 {
179 fprintf(debug, "{%d %d}", *i1+bFW ? dx : dy, *i2+bFW ? dy : dx );
180 }
181 if (bFW)
182 {
183 *i1 += dx;
184 *i2 += dy;
185 }
186 else
187 {
188 *i1 += dy;
189 *i2 += dx;
190 }
191 }
192
193 return cmp;
194}
195
196static int find_next_match_res(int *rnr1, int isize1,
197 int index1[], t_resinfo *resinfo1,
198 int *rnr2, int isize2,
199 int index2[], t_resinfo *resinfo2)
200{
201 int dx, dy, dmax, cmp, rr1, rr2;
202 gmx_bool bFW = FALSE0, bFF = FALSE0;
203
204 dx = dy = 0;
205 rr1 = 0;
206 while (rr1 < isize1 && *rnr1 != index1[rr1])
207 {
208 rr1++;
209 }
210 rr2 = 0;
211 while (rr2 < isize2 && *rnr2 != index2[rr2])
212 {
213 rr2++;
214 }
215
216 cmp = NOTSET-12345;
217 dmax = max(isize1-rr1, isize2-rr2)(((isize1-rr1) > (isize2-rr2)) ? (isize1-rr1) : (isize2-rr2
) );
218 if (debug)
219 {
220 fprintf(debug, " R:%d-%d:%d-%d:%d ",
221 rr1, isize1, rr2, isize2, dmax);
222 }
223 for (dx = 0; dx < dmax && cmp != 0; dx++)
224 {
225 for (dy = 0; dy <= dx && cmp != 0; dy++)
226 {
227 if (dx != dy)
228 {
229 cmp = NOTSET-12345;
230 if (rr1+dx < isize1 && rr2+dy < isize2)
231 {
232 bFW = TRUE1;
233 cmp = debug_strcmp(*resinfo1[index1[rr1+dx]].name,
234 *resinfo2[index2[rr2+dy]].name);
235 if (debug)
236 {
237 fprintf(debug, "(%d %d)", rr1+dx, rr2+dy);
238 }
239 }
240 if (cmp != 0 && rr1+dy < isize1 && rr2+dx < isize2)
241 {
242 bFW = FALSE0;
243 cmp = debug_strcmp(*resinfo1[index1[rr1+dy]].name,
244 *resinfo2[index2[rr2+dx]].name);
245 if (debug)
246 {
247 fprintf(debug, "(%d %d)", rr1+dy, rr2+dx);
248 }
249 }
250 if (dx != 0 && cmp != 0 && rr1+dx < isize1 && rr2+dx < isize2)
251 {
252 bFF = TRUE1;
253 cmp = debug_strcmp(*resinfo1[index1[rr1+dx]].name,
254 *resinfo2[index2[rr2+dx]].name);
255 if (debug)
256 {
257 fprintf(debug, "(%d %d)", rr1+dx, rr2+dx);
258 }
259 }
260 else
261 {
262 bFF = FALSE0;
263 }
264 }
265 }
266 }
267 /* apparently, dx and dy are incremented one more time
268 as the loops terminate: we correct this here */
269 dx--;
270 dy--;
271 /* if we skipped equal on both sides, only skip one residue
272 to allow for single mutations of residues... */
273 if (bFF)
274 {
275 if (debug)
276 {
277 fprintf(debug, "%d.%d.%dX%sX%s", dx, rr1, rr2,
278 *resinfo1[index1[rr1+1]].name,
279 *resinfo2[index2[rr2+1]].name);
280 }
281 dx = 1;
282 }
283 if (cmp == 0)
284 {
285 if (debug)
286 {
287 fprintf(debug, "!");
288 }
289 if (bFF)
290 {
291 rr1 += dx;
292 rr2 += dx;
293 }
294 else
295 if (bFW)
296 {
297 rr1 += dx;
298 rr2 += dy;
299 }
300 else
301 {
302 rr1 += dy;
303 rr2 += dx;
304 }
305 *rnr1 = index1[rr1];
306 *rnr2 = index2[rr2];
307 }
308
309 return cmp;
310}
311
312int find_first_atom_in_res(int rnr, int isize, atom_id index[], t_atom atom[])
313{
314 int i;
315
316 i = 0;
317 while (i < isize && atom[index[i]].resind != rnr)
318 {
319 i++;
320 }
321
322 if (atom[index[i]].resind == rnr)
323 {
324 return i;
325 }
326 else
327 {
328 return NOTSET-12345;
329 }
330}
331
332void find_matching_names(int *isize1, atom_id index1[], t_atoms *atoms1,
333 int *isize2, atom_id index2[], t_atoms *atoms2)
334{
335 int i, i1, i2, ii1, ii2, m1, m2;
336 int atcmp, rescmp;
337 int r, rnr1, rnr2, prnr1, prnr2;
338 int rsize1, rsize2;
339 int *rindex1, *rindex2;
340 char *resnm1, *resnm2, *atnm1, *atnm2;
341 char ***atnms1, ***atnms2;
342 t_resinfo *resinfo1, *resinfo2;
343
344 /* set some handy shortcuts */
345 resinfo1 = atoms1->resinfo;
346 atnms1 = atoms1->atomname;
347 resinfo2 = atoms2->resinfo;
348 atnms2 = atoms2->atomname;
349
350 /* build indexes of selected residues */
351 snew(rindex1, atoms1->nres)(rindex1) = save_calloc("rindex1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 351, (atoms1->nres), sizeof(*(rindex1)));
352 rsize1 = build_res_index(*isize1, index1, atoms1->atom, rindex1);
353 snew(rindex2, atoms2->nres)(rindex2) = save_calloc("rindex2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 353, (atoms2->nres), sizeof(*(rindex2)));
354 rsize2 = build_res_index(*isize2, index2, atoms2->atom, rindex2);
355
356 i1 = i2 = 0;
357 ii1 = ii2 = 0;
358 atcmp = rescmp = 0;
359 prnr1 = prnr2 = NOTSET-12345;
360 if (debug)
361 {
362 fprintf(debug, "Find matching names: %d, %d\n", *isize1, *isize2);
363 }
364 while (atcmp == 0 && i1 < *isize1 && i2 < *isize2)
365 {
366 /* prologue */
367 rnr1 = atoms1->atom[index1[i1]].resind;
368 rnr2 = atoms2->atom[index2[i2]].resind;
369 if (rnr1 != prnr1 || rnr2 != prnr2)
370 {
371 if (debug)
372 {
373 fprintf(debug, "R: %s%d %s%d\n",
374 *resinfo1[rnr1].name, rnr1, *resinfo2[rnr2].name, rnr2);
375 }
376 rescmp = strcmp(*resinfo1[rnr1].name, *resinfo2[rnr2].name)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(*resinfo1[rnr1].name) && __builtin_constant_p (*resinfo2
[rnr2].name) && (__s1_len = strlen (*resinfo1[rnr1].name
), __s2_len = strlen (*resinfo2[rnr2].name), (!((size_t)(const
void *)((*resinfo1[rnr1].name) + 1) - (size_t)(const void *)
(*resinfo1[rnr1].name) == 1) || __s1_len >= 4) && (
!((size_t)(const void *)((*resinfo2[rnr2].name) + 1) - (size_t
)(const void *)(*resinfo2[rnr2].name) == 1) || __s2_len >=
4)) ? __builtin_strcmp (*resinfo1[rnr1].name, *resinfo2[rnr2
].name) : (__builtin_constant_p (*resinfo1[rnr1].name) &&
((size_t)(const void *)((*resinfo1[rnr1].name) + 1) - (size_t
)(const void *)(*resinfo1[rnr1].name) == 1) && (__s1_len
= strlen (*resinfo1[rnr1].name), __s1_len < 4) ? (__builtin_constant_p
(*resinfo2[rnr2].name) && ((size_t)(const void *)((*
resinfo2[rnr2].name) + 1) - (size_t)(const void *)(*resinfo2[
rnr2].name) == 1) ? __builtin_strcmp (*resinfo1[rnr1].name, *
resinfo2[rnr2].name) : (__extension__ ({ const unsigned char *
__s2 = (const unsigned char *) (const char *) (*resinfo2[rnr2
].name); int __result = (((const unsigned char *) (const char
*) (*resinfo1[rnr1].name))[0] - __s2[0]); if (__s1_len > 0
&& __result == 0) { __result = (((const unsigned char
*) (const char *) (*resinfo1[rnr1].name))[1] - __s2[1]); if (
__s1_len > 1 && __result == 0) { __result = (((const
unsigned char *) (const char *) (*resinfo1[rnr1].name))[2] -
__s2[2]); if (__s1_len > 2 && __result == 0) __result
= (((const unsigned char *) (const char *) (*resinfo1[rnr1].
name))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
(*resinfo2[rnr2].name) && ((size_t)(const void *)((*
resinfo2[rnr2].name) + 1) - (size_t)(const void *)(*resinfo2[
rnr2].name) == 1) && (__s2_len = strlen (*resinfo2[rnr2
].name), __s2_len < 4) ? (__builtin_constant_p (*resinfo1[
rnr1].name) && ((size_t)(const void *)((*resinfo1[rnr1
].name) + 1) - (size_t)(const void *)(*resinfo1[rnr1].name) ==
1) ? __builtin_strcmp (*resinfo1[rnr1].name, *resinfo2[rnr2]
.name) : (- (__extension__ ({ const unsigned char *__s2 = (const
unsigned char *) (const char *) (*resinfo1[rnr1].name); int __result
= (((const unsigned char *) (const char *) (*resinfo2[rnr2].
name))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
*resinfo2[rnr2].name))[1] - __s2[1]); if (__s2_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (*resinfo2[rnr2].name))[2] - __s2[2]); if (__s2_len >
2 && __result == 0) __result = (((const unsigned char
*) (const char *) (*resinfo2[rnr2].name))[3] - __s2[3]); } }
__result; })))) : __builtin_strcmp (*resinfo1[rnr1].name, *resinfo2
[rnr2].name)))); });
377 }
378 if (debug)
379 {
380 fprintf(debug, "comparing %d %d", i1, i2);
381 }
382 atcmp = debug_strcmp(*atnms1[index1[i1]], *atnms2[index2[i2]]);
383
384 /* the works */
385 if (atcmp != 0) /* no match -> find match within residues */
386 {
387 m1 = find_res_end(i1, *isize1, index1, atoms1);
388 m2 = find_res_end(i2, *isize2, index2, atoms2);
389 if (debug)
390 {
391 fprintf(debug, " [%d<%d %d<%d]", i1, m1, i2, m2);
392 }
393 atcmp = find_next_match_atoms_in_res(&i1, index1, m1, atnms1,
394 &i2, index2, m2, atnms2);
395 if (debug)
396 {
397 fprintf(debug, " -> %d %d %s-%s", i1, i2,
398 *atnms1[index1[i1]], *atnms2[index2[i2]]);
399 }
400
401 }
402 if (atcmp != 0) /* still no match -> next residue(s) */
403 {
404 prnr1 = rnr1;
405 prnr2 = rnr2;
406 rescmp = find_next_match_res(&rnr1, rsize1, rindex1, resinfo1,
407 &rnr2, rsize2, rindex2, resinfo2);
408 if (rnr1 != prnr1)
409 {
410 i1 = find_first_atom_in_res(rnr1, *isize1, index1, atoms1->atom);
411 }
412 if (rnr2 != prnr2)
413 {
414 i2 = find_first_atom_in_res(rnr2, *isize2, index2, atoms2->atom);
415 }
416 if (debug)
417 {
418 fprintf(debug, " -> %s%d-%s%d %s%d-%s%d",
419 *resinfo1[rnr1].name, rnr1,
420 *atnms1[index1[i1]], index1[i1],
421 *resinfo2[rnr2].name, rnr2,
422 *atnms2[index2[i2]], index2[i2]);
423 }
424 m1 = find_res_end(i1, *isize1, index1, atoms1);
425 m2 = find_res_end(i2, *isize2, index2, atoms2);
426 if (debug)
427 {
428 fprintf(debug, " [%d<%d %d<%d]", i1, m1, i2, m2);
429 }
430 atcmp = find_next_match_atoms_in_res(&i1, index1, m1, atnms1,
431 &i2, index2, m2, atnms2);
432 if (debug)
433 {
434 fprintf(debug, " -> %d %d %s-%s", i1, i2,
435 *atnms1[index1[i1]], *atnms2[index2[i2]]);
436 }
437 }
438 if (debug)
439 {
440 fprintf(debug, "(%d %d): %d %d\n", ii1, ii2, atcmp, rescmp);
441 }
442 if (atcmp == 0) /* if match -> store indices */
443 {
444 index1[ii1++] = index1[i1];
445 index2[ii2++] = index2[i2];
446 }
447 i1++;
448 i2++;
449
450 /* epilogue */
451 prnr1 = rnr1;
452 prnr2 = rnr2;
453 }
454
455 if (ii1 != ii2)
456 {
457 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 457, "DEATH HORROR: non-equal number of matching atoms!\n");
458 }
459 if (ii1 == i1 && ii2 == i2)
460 {
461 printf("All atoms in index groups 1 and 2 match\n");
462 }
463 else
464 {
465 if (i1 == i2 && ii1 == ii2)
466 {
467 printf("Both index groups modified from %d to %d atoms\n", i1, ii1);
468 }
469 else
470 {
471 if (ii1 != i1)
472 {
473 printf("Index group 1 modified from %d to %d atoms\n", i1, ii1);
474 }
475 if (ii2 != i2)
476 {
477 printf("Index group 2 modified from %d to %d atoms\n", i2, ii2);
478 }
479 }
480 *isize1 = ii1;
481 *isize2 = ii2;
482 }
483}
484/* 1 */
485
486int gmx_confrms(int argc, char *argv[])
487{
488 const char *desc[] = {
489 "[THISMODULE] computes the root mean square deviation (RMSD) of two",
490 "structures after least-squares fitting the second structure on the first one.",
491 "The two structures do NOT need to have the same number of atoms,",
492 "only the two index groups used for the fit need to be identical.",
493 "With [TT]-name[tt] only matching atom names from the selected groups",
494 "will be used for the fit and RMSD calculation. This can be useful ",
495 "when comparing mutants of a protein.",
496 "[PAR]",
497 "The superimposed structures are written to file. In a [TT].pdb[tt] file",
498 "the two structures will be written as separate models",
499 "(use [TT]rasmol -nmrpdb[tt]). Also in a [TT].pdb[tt] file, B-factors",
500 "calculated from the atomic MSD values can be written with [TT]-bfac[tt].",
501 };
502 static gmx_bool bOne = FALSE0, bRmpbc = FALSE0, bMW = TRUE1, bName = FALSE0,
503 bBfac = FALSE0, bFit = TRUE1, bLabel = FALSE0;
504
505 t_pargs pa[] = {
506 { "-one", FALSE0, etBOOL, {&bOne}, "Only write the fitted structure to file" },
507 { "-mw", FALSE0, etBOOL, {&bMW}, "Mass-weighted fitting and RMSD" },
508 { "-pbc", FALSE0, etBOOL, {&bRmpbc}, "Try to make molecules whole again" },
509 { "-fit", FALSE0, etBOOL, {&bFit},
510 "Do least squares superposition of the target structure to the reference" },
511 { "-name", FALSE0, etBOOL, {&bName},
512 "Only compare matching atom names" },
513 { "-label", FALSE0, etBOOL, {&bLabel},
514 "Added chain labels A for first and B for second structure"},
515 { "-bfac", FALSE0, etBOOL, {&bBfac},
516 "Output B-factors from atomic MSD values" }
517 };
518 t_filenm fnm[] = {
519 { efTPS, "-f1", "conf1.gro", ffREAD1<<1 },
520 { efSTX, "-f2", "conf2", ffREAD1<<1 },
521 { efSTO, "-o", "fit.pdb", ffWRITE1<<2 },
522 { efNDX, "-n1", "fit1.ndx", ffOPTRD(1<<1 | 1<<3) },
523 { efNDX, "-n2", "fit2.ndx", ffOPTRD(1<<1 | 1<<3) },
524 { efNDX, "-no", "match.ndx", ffOPTWR(1<<2| 1<<3) }
525 };
526#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
527
528 /* the two structure files */
529 const char *conf1file, *conf2file, *matchndxfile, *outfile;
530 FILE *fp;
531 char title1[STRLEN4096], title2[STRLEN4096], *name1, *name2;
532 t_topology *top1, *top2;
533 int ePBC1, ePBC2;
534 t_atoms *atoms1, *atoms2;
535 int warn = 0;
536 atom_id at;
537 real *w_rls, mass, totmass;
538 rvec *x1, *v1, *x2, *v2, *fit_x;
539 matrix box1, box2;
540
541 output_env_t oenv;
542
543 /* counters */
544 int i, j, m;
545
546 /* center of mass calculation */
547 real tmas1, tmas2;
548 rvec xcm1, xcm2;
549
550 /* variables for fit */
551 char *groupnames1, *groupnames2;
552 int isize1, isize2;
553 atom_id *index1, *index2;
554 real rms, msd, minmsd, maxmsd;
555 real *msds;
556
557
558 if (!parse_common_args(&argc, argv, PCA_BE_NICE(1<<13) | PCA_CAN_VIEW(1<<5),
559 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))
560 {
561 return 0;
562 }
563 matchndxfile = opt2fn_null("-no", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
564 conf1file = ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
565 conf2file = ftp2fn(efSTX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
566
567 /* reading reference structure from first structure file */
568 fprintf(stderrstderr, "\nReading first structure file\n");
569 snew(top1, 1)(top1) = save_calloc("top1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 569, (1), sizeof(*(top1)));
570 read_tps_conf(conf1file, title1, top1, &ePBC1, &x1, &v1, box1, TRUE1);
571 atoms1 = &(top1->atoms);
572 fprintf(stderrstderr, "%s\nContaining %d atoms in %d residues\n",
573 title1, atoms1->nr, atoms1->nres);
574
575 if (bRmpbc)
576 {
577 rm_gropbc(atoms1, x1, box1);
578 }
579
580 fprintf(stderrstderr, "Select group from first structure\n");
581 get_index(atoms1, opt2fn_null("-n1", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
582 1, &isize1, &index1, &groupnames1);
583 printf("\n");
584
585 if (bFit && (isize1 < 3))
586 {
587 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 587, "Need >= 3 points to fit!\n");
588 }
589
590 /* reading second structure file */
591 fprintf(stderrstderr, "\nReading second structure file\n");
592 snew(top2, 1)(top2) = save_calloc("top2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 592, (1), sizeof(*(top2)));
593 read_tps_conf(conf2file, title2, top2, &ePBC2, &x2, &v2, box2, TRUE1);
594 atoms2 = &(top2->atoms);
595 fprintf(stderrstderr, "%s\nContaining %d atoms in %d residues\n",
596 title2, atoms2->nr, atoms2->nres);
597
598 if (bRmpbc)
599 {
600 rm_gropbc(atoms2, x2, box2);
601 }
602
603 fprintf(stderrstderr, "Select group from second structure\n");
604 get_index(atoms2, opt2fn_null("-n2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
605 1, &isize2, &index2, &groupnames2);
606
607 if (bName)
608 {
609 find_matching_names(&isize1, index1, atoms1, &isize2, index2, atoms2);
610 if (matchndxfile)
611 {
612 fp = gmx_ffopen(matchndxfile, "w");
613 fprintf(fp, "; Matching atoms between %s from %s and %s from %s\n",
614 groupnames1, conf1file, groupnames2, conf2file);
615 fprintf(fp, "[ Match_%s_%s ]\n", conf1file, groupnames1);
616 for (i = 0; i < isize1; i++)
617 {
618 fprintf(fp, "%4u%s", index1[i]+1, (i%15 == 14 || i == isize1-1) ? "\n" : " ");
619 }
620 fprintf(fp, "[ Match_%s_%s ]\n", conf2file, groupnames2);
621 for (i = 0; i < isize2; i++)
622 {
623 fprintf(fp, "%4u%s", index2[i]+1, (i%15 == 14 || i == isize2-1) ? "\n" : " ");
624 }
625 }
626 }
627
628 /* check isizes, must be equal */
629 if (isize2 != isize1)
630 {
631 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 631, "You selected groups with differen number of atoms.\n");
632 }
633
634 for (i = 0; i < isize1; i++)
635 {
636 name1 = *atoms1->atomname[index1[i]];
637 name2 = *atoms2->atomname[index2[i]];
638 if (strcmp(name1, name2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(name1) && __builtin_constant_p (name2) && (
__s1_len = strlen (name1), __s2_len = strlen (name2), (!((size_t
)(const void *)((name1) + 1) - (size_t)(const void *)(name1) ==
1) || __s1_len >= 4) && (!((size_t)(const void *)
((name2) + 1) - (size_t)(const void *)(name2) == 1) || __s2_len
>= 4)) ? __builtin_strcmp (name1, name2) : (__builtin_constant_p
(name1) && ((size_t)(const void *)((name1) + 1) - (size_t
)(const void *)(name1) == 1) && (__s1_len = strlen (name1
), __s1_len < 4) ? (__builtin_constant_p (name2) &&
((size_t)(const void *)((name2) + 1) - (size_t)(const void *
)(name2) == 1) ? __builtin_strcmp (name1, name2) : (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (name2); int __result = (((const unsigned char *) (const
char *) (name1))[0] - __s2[0]); if (__s1_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (name1))[1] - __s2[1]); if (__s1_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (name1))[2] - __s2[2]); if (__s1_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (name1))[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p
(name2) && ((size_t)(const void *)((name2) + 1) - (size_t
)(const void *)(name2) == 1) && (__s2_len = strlen (name2
), __s2_len < 4) ? (__builtin_constant_p (name1) &&
((size_t)(const void *)((name1) + 1) - (size_t)(const void *
)(name1) == 1) ? __builtin_strcmp (name1, name2) : (- (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (name1); int __result = (((const unsigned char *) (const
char *) (name2))[0] - __s2[0]); if (__s2_len > 0 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (name2))[1] - __s2[1]); if (__s2_len > 1 &&
__result == 0) { __result = (((const unsigned char *) (const
char *) (name2))[2] - __s2[2]); if (__s2_len > 2 &&
__result == 0) __result = (((const unsigned char *) (const char
*) (name2))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp
(name1, name2)))); }))
639 {
640 if (warn < 20)
641 {
642 fprintf(stderrstderr,
643 "Warning: atomnames at index %d don't match: %u %s, %u %s\n",
644 i+1, index1[i]+1, name1, index2[i]+1, name2);
645 }
646 warn++;
647 }
648 if (!bMW)
649 {
650 atoms1->atom[index1[i]].m = 1;
651 atoms2->atom[index2[i]].m = 1;
652 }
653 }
654 if (warn)
655 {
656 fprintf(stderrstderr, "%d atomname%s did not match\n", warn, (warn == 1) ? "" : "s");
657 }
658
659 if (bFit)
660 {
661 /* calculate and remove center of mass of structures */
662 calc_rm_cm(isize1, index1, atoms1, x1, xcm1);
663 calc_rm_cm(isize2, index2, atoms2, x2, xcm2);
664
665 snew(w_rls, atoms2->nr)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 665, (atoms2->nr), sizeof(*(w_rls)));
666 snew(fit_x, atoms2->nr)(fit_x) = save_calloc("fit_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 666, (atoms2->nr), sizeof(*(fit_x)));
667 for (at = 0; (at < isize1); at++)
668 {
669 w_rls[index2[at]] = atoms1->atom[index1[at]].m;
670 copy_rvec(x1[index1[at]], fit_x[index2[at]]);
671 }
672
673 /* do the least squares fit to the reference structure */
674 do_fit(atoms2->nr, w_rls, fit_x, x2);
675
676 sfree(fit_x)save_free("fit_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 676, (fit_x));
677 sfree(w_rls)save_free("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 677, (w_rls));
678 w_rls = NULL((void*)0);
679 }
680 else
681 {
682 clear_rvec(xcm1);
683 clear_rvec(xcm2);
684 w_rls = NULL((void*)0);
685 }
686
687 /* calculate the rms deviation */
688 rms = 0;
689 totmass = 0;
690 maxmsd = -1e18;
691 minmsd = 1e18;
692 snew(msds, isize1)(msds) = save_calloc("msds", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 692, (isize1), sizeof(*(msds)));
693 for (at = 0; at < isize1; at++)
694 {
695 mass = atoms1->atom[index1[at]].m;
696 for (m = 0; m < DIM3; m++)
697 {
698 msd = sqr(x1[index1[at]][m] - x2[index2[at]][m]);
699 rms += msd*mass;
700 msds[at] += msd;
701 }
702 maxmsd = max(maxmsd, msds[at])(((maxmsd) > (msds[at])) ? (maxmsd) : (msds[at]) );
703 minmsd = min(minmsd, msds[at])(((minmsd) < (msds[at])) ? (minmsd) : (msds[at]) );
704 totmass += mass;
705 }
706 rms = sqrt(rms/totmass);
707
708 printf("Root mean square deviation after lsq fit = %g nm\n", rms);
709 if (bBfac)
710 {
711 printf("Atomic MSD's range from %g to %g nm^2\n", minmsd, maxmsd);
712 }
713
714 if (bFit)
715 {
716 /* reset coordinates of reference and fitted structure */
717 for (i = 0; i < atoms1->nr; i++)
718 {
719 for (m = 0; m < DIM3; m++)
720 {
721 x1[i][m] += xcm1[m];
722 }
723 }
724 for (i = 0; i < atoms2->nr; i++)
725 {
726 for (m = 0; m < DIM3; m++)
727 {
728 x2[i][m] += xcm1[m];
729 }
730 }
731 }
732
733 outfile = ftp2fn(efSTO, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
734 switch (fn2ftp(outfile))
735 {
736 case efPDB:
737 case efBRK:
738 case efENT:
739 if (bBfac || bLabel)
740 {
741 srenew(atoms1->pdbinfo, atoms1->nr)(atoms1->pdbinfo) = save_realloc("atoms1->pdbinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 741, (atoms1->pdbinfo), (atoms1->nr), sizeof(*(atoms1
->pdbinfo)));
742 srenew(atoms1->atom, atoms1->nr)(atoms1->atom) = save_realloc("atoms1->atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 742, (atoms1->atom), (atoms1->nr), sizeof(*(atoms1->
atom))); /* Why renew atom? */
743
744 /* Avoid segfaults when writing the pdb-file */
745 for (i = 0; i < atoms1->nr; i++)
746 {
747 atoms1->pdbinfo[i].type = eptAtom;
748 atoms1->pdbinfo[i].occup = 1.00;
749 atoms1->pdbinfo[i].bAnisotropic = FALSE0;
750 if (bBfac)
751 {
752 atoms1->pdbinfo[i].bfac = 0;
753 }
754 if (bLabel)
755 {
756 atoms1->resinfo[atoms1->atom[i].resind].chainid = 'A';
757 }
758 }
759
760 for (i = 0; i < isize1; i++)
761 {
762 /* atoms1->pdbinfo[index1[i]].type = eptAtom; */
763/* atoms1->pdbinfo[index1[i]].bAnisotropic = FALSE; */
764 if (bBfac)
765 {
766 atoms1->pdbinfo[index1[i]].bfac = (800*M_PI3.14159265358979323846*M_PI3.14159265358979323846/3.0)*msds[i];
767 }
768/* if (bLabel) */
769/* atoms1->resinfo[atoms1->atom[index1[i]].resind].chain = 'A'; */
770 }
771 srenew(atoms2->pdbinfo, atoms2->nr)(atoms2->pdbinfo) = save_realloc("atoms2->pdbinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 771, (atoms2->pdbinfo), (atoms2->nr), sizeof(*(atoms2
->pdbinfo)));
772 srenew(atoms2->atom, atoms2->nr)(atoms2->atom) = save_realloc("atoms2->atom", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_confrms.c"
, 772, (atoms2->atom), (atoms2->nr), sizeof(*(atoms2->
atom))); /* Why renew atom? */
773
774 for (i = 0; i < atoms2->nr; i++)
775 {
776 atoms2->pdbinfo[i].type = eptAtom;
777 atoms2->pdbinfo[i].occup = 1.00;
778 atoms2->pdbinfo[i].bAnisotropic = FALSE0;
779 if (bBfac)
780 {
781 atoms2->pdbinfo[i].bfac = 0;
782 }
783 if (bLabel)
784 {
785 atoms2->resinfo[atoms1->atom[i].resind].chainid = 'B';
786 }
787 }
788
789 for (i = 0; i < isize2; i++)
790 {
791 /* atoms2->pdbinfo[index2[i]].type = eptAtom; */
792/* atoms2->pdbinfo[index2[i]].bAnisotropic = FALSE; */
793 if (bBfac)
794 {
795 atoms2->pdbinfo[index2[i]].bfac = (800*M_PI3.14159265358979323846*M_PI3.14159265358979323846/3.0)*msds[i];
796 }
797/* if (bLabel) */
798/* atoms2->resinfo[atoms2->atom[index2[i]].resind].chain = 'B'; */
799 }
800 }
801 fp = gmx_ffopen(outfile, "w");
802 if (!bOne)
803 {
804 write_pdbfile(fp, title1, atoms1, x1, ePBC1, box1, ' ', 1, NULL((void*)0), TRUE1);
805 }
806 write_pdbfile(fp, title2, atoms2, x2, ePBC2, box2, ' ', bOne ? -1 : 2, NULL((void*)0), TRUE1);
807 gmx_ffclose(fp);
808 break;
809 case efGRO:
810 if (bBfac)
811 {
812 fprintf(stderrstderr, "WARNING: cannot write B-factor values to gro file\n");
813 }
814 fp = gmx_ffopen(outfile, "w");
815 if (!bOne)
816 {
817 write_hconf_p(fp, title1, atoms1, 3, x1, v1, box1);
818 }
819 write_hconf_p(fp, title2, atoms2, 3, x2, v2, box2);
820 gmx_ffclose(fp);
821 break;
822 default:
823 if (bBfac)
824 {
825 fprintf(stderrstderr, "WARNING: cannot write B-factor values to %s file\n",
826 ftp2ext(fn2ftp(outfile)));
827 }
828 if (!bOne)
829 {
830 fprintf(stderrstderr,
831 "WARNING: cannot write the reference structure to %s file\n",
832 ftp2ext(fn2ftp(outfile)));
833 }
834 write_sto_conf(outfile, title2, atoms2, x2, v2, ePBC2, box2);
835 break;
836 }
837
838 view_all(oenv, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
839
840 return 0;
841}