Bug Summary

File:gromacs/gmxana/gmx_rmsdist.c
Location:line 790, column 5
Description:Value stored to 'natom' 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
43#include "macros.h"
44#include "gromacs/utility/smalloc.h"
45#include "typedefs.h"
46#include "names.h"
47#include "gromacs/commandline/pargs.h"
48#include "gromacs/fileio/tpxio.h"
49#include "gromacs/fileio/trxio.h"
50#include "gromacs/fileio/strdb.h"
51#include "gromacs/math/vec.h"
52#include "macros.h"
53#include "index.h"
54#include "pbc.h"
55#include "gromacs/fileio/xvgr.h"
56#include "viewit.h"
57#include "gromacs/utility/futil.h"
58#include "gromacs/fileio/matio.h"
59#include "gmx_ana.h"
60
61
62static void calc_dist(int nind, atom_id index[], rvec x[], int ePBC, matrix box,
63 real **d)
64{
65 int i, j;
66 real *xi;
67 rvec dx;
68 real temp2;
69 t_pbc pbc;
70
71 set_pbc(&pbc, ePBC, box);
72 for (i = 0; (i < nind-1); i++)
73 {
74 xi = x[index[i]];
75 for (j = i+1; (j < nind); j++)
76 {
77 pbc_dx(&pbc, xi, x[index[j]], dx);
78 temp2 = norm2(dx);
79 d[i][j] = sqrt(temp2);
80 }
81 }
82}
83
84static void calc_dist_tot(int nind, atom_id index[], rvec x[],
85 int ePBC, matrix box,
86 real **d, real **dtot, real **dtot2,
87 gmx_bool bNMR, real **dtot1_3, real **dtot1_6)
88{
89 int i, j;
90 real *xi;
91 real temp, temp2, temp1_3;
92 rvec dx;
93 t_pbc pbc;
94
95 set_pbc(&pbc, ePBC, box);
96 for (i = 0; (i < nind-1); i++)
97 {
98 xi = x[index[i]];
99 for (j = i+1; (j < nind); j++)
100 {
101 pbc_dx(&pbc, xi, x[index[j]], dx);
102 temp2 = norm2(dx);
103 temp = sqrt(temp2);
104 d[i][j] = temp;
105 dtot[i][j] += temp;
106 dtot2[i][j] += temp2;
107 if (bNMR)
108 {
109 temp1_3 = 1.0/(temp*temp2);
110 dtot1_3[i][j] += temp1_3;
111 dtot1_6[i][j] += temp1_3*temp1_3;
112 }
113 }
114 }
115}
116
117static void calc_nmr(int nind, int nframes, real **dtot1_3, real **dtot1_6,
118 real *max1_3, real *max1_6)
119{
120 int i, j;
121 real temp1_3, temp1_6;
122
123 for (i = 0; (i < nind-1); i++)
124 {
125 for (j = i+1; (j < nind); j++)
126 {
127 temp1_3 = pow(dtot1_3[i][j]/nframes, -1.0/3.0);
128 temp1_6 = pow(dtot1_6[i][j]/nframes, -1.0/6.0);
129 if (temp1_3 > *max1_3)
130 {
131 *max1_3 = temp1_3;
132 }
133 if (temp1_6 > *max1_6)
134 {
135 *max1_6 = temp1_6;
136 }
137 dtot1_3[i][j] = temp1_3;
138 dtot1_6[i][j] = temp1_6;
139 dtot1_3[j][i] = temp1_3;
140 dtot1_6[j][i] = temp1_6;
141 }
142 }
143}
144
145static char Hnum[] = "123";
146
147typedef struct {
148 int nr;
149 real r_3;
150 real r_6;
151 real i_3;
152 real i_6;
153} t_noe;
154
155typedef struct {
156 int anr;
157 int ianr;
158 int rnr;
159 char *aname;
160 char *rname;
161} t_noe_gr;
162
163typedef struct {
164 int rnr;
165 char *nname;
166 char *rname;
167 char *aname;
168} t_equiv;
169
170static int read_equiv(const char *eq_fn, t_equiv ***equivptr)
171{
172 FILE *fp;
173 char line[STRLEN4096], resname[10], atomname[10], *lp;
174 int neq, na, n, resnr;
175 t_equiv **equiv;
176
177 fp = gmx_ffopen(eq_fn, "r");
178 neq = 0;
179 equiv = NULL((void*)0);
180 while (get_a_line(fp, line, STRLEN4096))
181 {
182 lp = line;
183 /* this is not efficient, but I'm lazy */
184 srenew(equiv, neq+1)(equiv) = save_realloc("equiv", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 184, (equiv), (neq+1), sizeof(*(equiv)));
185 equiv[neq] = NULL((void*)0);
186 na = 0;
187 if (sscanf(lp, "%s %n", atomname, &n) == 1)
188 {
189 lp += n;
190 snew(equiv[neq], 1)(equiv[neq]) = save_calloc("equiv[neq]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 190, (1), sizeof(*(equiv[neq])));
191 equiv[neq][0].nname = strdup(atomname)(__extension__ (__builtin_constant_p (atomname) && ((
size_t)(const void *)((atomname) + 1) - (size_t)(const void *
)(atomname) == 1) ? (((const char *) (atomname))[0] == '\0' ?
(char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len =
strlen (atomname) + 1; char *__retval = (char *) malloc (__len
); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval
, atomname, __len); __retval; })) : __strdup (atomname)));
192 while (sscanf(lp, "%d %s %s %n", &resnr, resname, atomname, &n) == 3)
193 {
194 /* this is not efficient, but I'm lazy (again) */
195 srenew(equiv[neq], na+1)(equiv[neq]) = save_realloc("equiv[neq]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 195, (equiv[neq]), (na+1), sizeof(*(equiv[neq])));
196 equiv[neq][na].rnr = resnr-1;
197 equiv[neq][na].rname = strdup(resname)(__extension__ (__builtin_constant_p (resname) && ((size_t
)(const void *)((resname) + 1) - (size_t)(const void *)(resname
) == 1) ? (((const char *) (resname))[0] == '\0' ? (char *) calloc
((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (resname
) + 1; char *__retval = (char *) malloc (__len); if (__retval
!= ((void*)0)) __retval = (char *) memcpy (__retval, resname
, __len); __retval; })) : __strdup (resname)));
198 equiv[neq][na].aname = strdup(atomname)(__extension__ (__builtin_constant_p (atomname) && ((
size_t)(const void *)((atomname) + 1) - (size_t)(const void *
)(atomname) == 1) ? (((const char *) (atomname))[0] == '\0' ?
(char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len =
strlen (atomname) + 1; char *__retval = (char *) malloc (__len
); if (__retval != ((void*)0)) __retval = (char *) memcpy (__retval
, atomname, __len); __retval; })) : __strdup (atomname)));
199 if (na > 0)
200 {
201 equiv[neq][na].nname = NULL((void*)0);
202 }
203 na++;
204 lp += n;
205 }
206 }
207 /* make empty element as flag for end of array */
208 srenew(equiv[neq], na+1)(equiv[neq]) = save_realloc("equiv[neq]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 208, (equiv[neq]), (na+1), sizeof(*(equiv[neq])));
209 equiv[neq][na].rnr = NOTSET-12345;
210 equiv[neq][na].rname = NULL((void*)0);
211 equiv[neq][na].aname = NULL((void*)0);
212
213 /* next */
214 neq++;
215 }
216 gmx_ffclose(fp);
217
218 *equivptr = equiv;
219
220 return neq;
221}
222
223static void dump_equiv(FILE *out, int neq, t_equiv **equiv)
224{
225 int i, j;
226
227 fprintf(out, "Dumping equivalent list\n");
228 for (i = 0; i < neq; i++)
229 {
230 fprintf(out, "%s", equiv[i][0].nname);
231 for (j = 0; equiv[i][j].rnr != NOTSET-12345; j++)
232 {
233 fprintf(out, " %d %s %s",
234 equiv[i][j].rnr, equiv[i][j].rname, equiv[i][j].aname);
235 }
236 fprintf(out, "\n");
237 }
238}
239
240static gmx_bool is_equiv(int neq, t_equiv **equiv, char **nname,
241 int rnr1, char *rname1, char *aname1,
242 int rnr2, char *rname2, char *aname2)
243{
244 int i, j;
245 gmx_bool bFound;
246
247 bFound = FALSE0;
248 /* we can terminate each loop when bFound is true! */
249 for (i = 0; i < neq && !bFound; i++)
250 {
251 /* find first atom */
252 for (j = 0; equiv[i][j].rnr != NOTSET-12345 && !bFound; j++)
253 {
254 bFound = ( equiv[i][j].rnr == rnr1 &&
255 strcmp(equiv[i][j].rname, rname1)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(equiv[i][j].rname) && __builtin_constant_p (rname1)
&& (__s1_len = strlen (equiv[i][j].rname), __s2_len =
strlen (rname1), (!((size_t)(const void *)((equiv[i][j].rname
) + 1) - (size_t)(const void *)(equiv[i][j].rname) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)((rname1) + 1) -
(size_t)(const void *)(rname1) == 1) || __s2_len >= 4)) ?
__builtin_strcmp (equiv[i][j].rname, rname1) : (__builtin_constant_p
(equiv[i][j].rname) && ((size_t)(const void *)((equiv
[i][j].rname) + 1) - (size_t)(const void *)(equiv[i][j].rname
) == 1) && (__s1_len = strlen (equiv[i][j].rname), __s1_len
< 4) ? (__builtin_constant_p (rname1) && ((size_t
)(const void *)((rname1) + 1) - (size_t)(const void *)(rname1
) == 1) ? __builtin_strcmp (equiv[i][j].rname, rname1) : (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (rname1); int __result = (((const unsigned char *) (
const char *) (equiv[i][j].rname))[0] - __s2[0]); if (__s1_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (equiv[i][j].rname))[1] - __s2[1]); if
(__s1_len > 1 && __result == 0) { __result = (((const
unsigned char *) (const char *) (equiv[i][j].rname))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (equiv[i][j].rname)
)[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (
rname1) && ((size_t)(const void *)((rname1) + 1) - (size_t
)(const void *)(rname1) == 1) && (__s2_len = strlen (
rname1), __s2_len < 4) ? (__builtin_constant_p (equiv[i][j
].rname) && ((size_t)(const void *)((equiv[i][j].rname
) + 1) - (size_t)(const void *)(equiv[i][j].rname) == 1) ? __builtin_strcmp
(equiv[i][j].rname, rname1) : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (equiv[i
][j].rname); int __result = (((const unsigned char *) (const char
*) (rname1))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
rname1))[1] - __s2[1]); if (__s2_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
rname1))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (rname1
))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (equiv
[i][j].rname, rname1)))); }) == 0 &&
256 strcmp(equiv[i][j].aname, aname1)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(equiv[i][j].aname) && __builtin_constant_p (aname1)
&& (__s1_len = strlen (equiv[i][j].aname), __s2_len =
strlen (aname1), (!((size_t)(const void *)((equiv[i][j].aname
) + 1) - (size_t)(const void *)(equiv[i][j].aname) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)((aname1) + 1) -
(size_t)(const void *)(aname1) == 1) || __s2_len >= 4)) ?
__builtin_strcmp (equiv[i][j].aname, aname1) : (__builtin_constant_p
(equiv[i][j].aname) && ((size_t)(const void *)((equiv
[i][j].aname) + 1) - (size_t)(const void *)(equiv[i][j].aname
) == 1) && (__s1_len = strlen (equiv[i][j].aname), __s1_len
< 4) ? (__builtin_constant_p (aname1) && ((size_t
)(const void *)((aname1) + 1) - (size_t)(const void *)(aname1
) == 1) ? __builtin_strcmp (equiv[i][j].aname, aname1) : (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (aname1); int __result = (((const unsigned char *) (
const char *) (equiv[i][j].aname))[0] - __s2[0]); if (__s1_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (equiv[i][j].aname))[1] - __s2[1]); if
(__s1_len > 1 && __result == 0) { __result = (((const
unsigned char *) (const char *) (equiv[i][j].aname))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (equiv[i][j].aname)
)[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (
aname1) && ((size_t)(const void *)((aname1) + 1) - (size_t
)(const void *)(aname1) == 1) && (__s2_len = strlen (
aname1), __s2_len < 4) ? (__builtin_constant_p (equiv[i][j
].aname) && ((size_t)(const void *)((equiv[i][j].aname
) + 1) - (size_t)(const void *)(equiv[i][j].aname) == 1) ? __builtin_strcmp
(equiv[i][j].aname, aname1) : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (equiv[i
][j].aname); int __result = (((const unsigned char *) (const char
*) (aname1))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
aname1))[1] - __s2[1]); if (__s2_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
aname1))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (aname1
))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (equiv
[i][j].aname, aname1)))); }) == 0 );
257 }
258 if (bFound)
259 {
260 /* find second atom */
261 bFound = FALSE0;
262 for (j = 0; equiv[i][j].rnr != NOTSET-12345 && !bFound; j++)
263 {
264 bFound = ( equiv[i][j].rnr == rnr2 &&
265 strcmp(equiv[i][j].rname, rname2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(equiv[i][j].rname) && __builtin_constant_p (rname2)
&& (__s1_len = strlen (equiv[i][j].rname), __s2_len =
strlen (rname2), (!((size_t)(const void *)((equiv[i][j].rname
) + 1) - (size_t)(const void *)(equiv[i][j].rname) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)((rname2) + 1) -
(size_t)(const void *)(rname2) == 1) || __s2_len >= 4)) ?
__builtin_strcmp (equiv[i][j].rname, rname2) : (__builtin_constant_p
(equiv[i][j].rname) && ((size_t)(const void *)((equiv
[i][j].rname) + 1) - (size_t)(const void *)(equiv[i][j].rname
) == 1) && (__s1_len = strlen (equiv[i][j].rname), __s1_len
< 4) ? (__builtin_constant_p (rname2) && ((size_t
)(const void *)((rname2) + 1) - (size_t)(const void *)(rname2
) == 1) ? __builtin_strcmp (equiv[i][j].rname, rname2) : (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (rname2); int __result = (((const unsigned char *) (
const char *) (equiv[i][j].rname))[0] - __s2[0]); if (__s1_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (equiv[i][j].rname))[1] - __s2[1]); if
(__s1_len > 1 && __result == 0) { __result = (((const
unsigned char *) (const char *) (equiv[i][j].rname))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (equiv[i][j].rname)
)[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (
rname2) && ((size_t)(const void *)((rname2) + 1) - (size_t
)(const void *)(rname2) == 1) && (__s2_len = strlen (
rname2), __s2_len < 4) ? (__builtin_constant_p (equiv[i][j
].rname) && ((size_t)(const void *)((equiv[i][j].rname
) + 1) - (size_t)(const void *)(equiv[i][j].rname) == 1) ? __builtin_strcmp
(equiv[i][j].rname, rname2) : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (equiv[i
][j].rname); int __result = (((const unsigned char *) (const char
*) (rname2))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
rname2))[1] - __s2[1]); if (__s2_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
rname2))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (rname2
))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (equiv
[i][j].rname, rname2)))); }) == 0 &&
266 strcmp(equiv[i][j].aname, aname2)__extension__ ({ size_t __s1_len, __s2_len; (__builtin_constant_p
(equiv[i][j].aname) && __builtin_constant_p (aname2)
&& (__s1_len = strlen (equiv[i][j].aname), __s2_len =
strlen (aname2), (!((size_t)(const void *)((equiv[i][j].aname
) + 1) - (size_t)(const void *)(equiv[i][j].aname) == 1) || __s1_len
>= 4) && (!((size_t)(const void *)((aname2) + 1) -
(size_t)(const void *)(aname2) == 1) || __s2_len >= 4)) ?
__builtin_strcmp (equiv[i][j].aname, aname2) : (__builtin_constant_p
(equiv[i][j].aname) && ((size_t)(const void *)((equiv
[i][j].aname) + 1) - (size_t)(const void *)(equiv[i][j].aname
) == 1) && (__s1_len = strlen (equiv[i][j].aname), __s1_len
< 4) ? (__builtin_constant_p (aname2) && ((size_t
)(const void *)((aname2) + 1) - (size_t)(const void *)(aname2
) == 1) ? __builtin_strcmp (equiv[i][j].aname, aname2) : (__extension__
({ const unsigned char *__s2 = (const unsigned char *) (const
char *) (aname2); int __result = (((const unsigned char *) (
const char *) (equiv[i][j].aname))[0] - __s2[0]); if (__s1_len
> 0 && __result == 0) { __result = (((const unsigned
char *) (const char *) (equiv[i][j].aname))[1] - __s2[1]); if
(__s1_len > 1 && __result == 0) { __result = (((const
unsigned char *) (const char *) (equiv[i][j].aname))[2] - __s2
[2]); if (__s1_len > 2 && __result == 0) __result =
(((const unsigned char *) (const char *) (equiv[i][j].aname)
)[3] - __s2[3]); } } __result; }))) : (__builtin_constant_p (
aname2) && ((size_t)(const void *)((aname2) + 1) - (size_t
)(const void *)(aname2) == 1) && (__s2_len = strlen (
aname2), __s2_len < 4) ? (__builtin_constant_p (equiv[i][j
].aname) && ((size_t)(const void *)((equiv[i][j].aname
) + 1) - (size_t)(const void *)(equiv[i][j].aname) == 1) ? __builtin_strcmp
(equiv[i][j].aname, aname2) : (- (__extension__ ({ const unsigned
char *__s2 = (const unsigned char *) (const char *) (equiv[i
][j].aname); int __result = (((const unsigned char *) (const char
*) (aname2))[0] - __s2[0]); if (__s2_len > 0 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
aname2))[1] - __s2[1]); if (__s2_len > 1 && __result
== 0) { __result = (((const unsigned char *) (const char *) (
aname2))[2] - __s2[2]); if (__s2_len > 2 && __result
== 0) __result = (((const unsigned char *) (const char *) (aname2
))[3] - __s2[3]); } } __result; })))) : __builtin_strcmp (equiv
[i][j].aname, aname2)))); }) == 0 );
267 }
268 }
269 }
270 if (bFound)
271 {
272 *nname = strdup(equiv[i-1][0].nname)(__extension__ (__builtin_constant_p (equiv[i-1][0].nname) &&
((size_t)(const void *)((equiv[i-1][0].nname) + 1) - (size_t
)(const void *)(equiv[i-1][0].nname) == 1) ? (((const char *)
(equiv[i-1][0].nname))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (equiv[i-1][0].nname
) + 1; char *__retval = (char *) malloc (__len); if (__retval
!= ((void*)0)) __retval = (char *) memcpy (__retval, equiv[i
-1][0].nname, __len); __retval; })) : __strdup (equiv[i-1][0]
.nname)));
273 }
274
275 return bFound;
276}
277
278static int analyze_noe_equivalent(const char *eq_fn,
279 t_atoms *atoms, int isize, atom_id *index,
280 gmx_bool bSumH,
281 atom_id *noe_index, t_noe_gr *noe_gr)
282{
283 FILE *fp;
284 int i, j, anmil, anmjl, rnri, rnrj, gi, groupnr, neq;
285 char *anmi, *anmj, **nnm;
286 gmx_bool bMatch, bEquiv;
287 t_equiv **equiv;
288
289 snew(nnm, isize)(nnm) = save_calloc("nnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 289, (isize), sizeof(*(nnm)));
290 if (bSumH)
291 {
292 if (eq_fn)
293 {
294 neq = read_equiv(eq_fn, &equiv);
295 if (debug)
296 {
297 dump_equiv(debug, neq, equiv);
298 }
299 }
300 else
301 {
302 neq = 0;
303 equiv = NULL((void*)0);
304 }
305
306 groupnr = 0;
307 for (i = 0; i < isize; i++)
308 {
309 if (equiv && i < isize-1)
310 {
311 /* check explicit list of equivalent atoms */
312 do
313 {
314 j = i+1;
315 rnri = atoms->atom[index[i]].resind;
316 rnrj = atoms->atom[index[j]].resind;
317 bEquiv =
318 is_equiv(neq, equiv, &nnm[i],
319 rnri, *atoms->resinfo[rnri].name, *atoms->atomname[index[i]],
320 rnrj, *atoms->resinfo[rnrj].name, *atoms->atomname[index[j]]);
321 if (nnm[i] && bEquiv)
322 {
323 nnm[j] = strdup(nnm[i])(__extension__ (__builtin_constant_p (nnm[i]) && ((size_t
)(const void *)((nnm[i]) + 1) - (size_t)(const void *)(nnm[i]
) == 1) ? (((const char *) (nnm[i]))[0] == '\0' ? (char *) calloc
((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (nnm[i])
+ 1; char *__retval = (char *) malloc (__len); if (__retval !=
((void*)0)) __retval = (char *) memcpy (__retval, nnm[i], __len
); __retval; })) : __strdup (nnm[i])));
324 }
325 if (bEquiv)
326 {
327 /* set index for matching atom */
328 noe_index[j] = groupnr;
329 /* skip matching atom */
330 i = j;
331 }
332 }
333 while (bEquiv && i < isize-1);
334 }
335 else
336 {
337 bEquiv = FALSE0;
338 }
339 if (!bEquiv)
340 {
341 /* look for triplets of consecutive atoms with name XX?,
342 X are any number of letters or digits and ? goes from 1 to 3
343 This is supposed to cover all CH3 groups and the like */
344 anmi = *atoms->atomname[index[i]];
345 anmil = strlen(anmi);
346 bMatch = i < isize-3 && anmi[anmil-1] == '1';
347 if (bMatch)
348 {
349 for (j = 1; j < 3; j++)
350 {
351 anmj = *atoms->atomname[index[i+j]];
352 anmjl = strlen(anmj);
353 bMatch = bMatch && ( anmil == anmjl && anmj[anmjl-1] == Hnum[j] &&
354 strncmp(anmi, anmj, anmil-1)(__extension__ (__builtin_constant_p (anmil-1) && ((__builtin_constant_p
(anmi) && strlen (anmi) < ((size_t) (anmil-1))) ||
(__builtin_constant_p (anmj) && strlen (anmj) < (
(size_t) (anmil-1)))) ? __extension__ ({ size_t __s1_len, __s2_len
; (__builtin_constant_p (anmi) && __builtin_constant_p
(anmj) && (__s1_len = strlen (anmi), __s2_len = strlen
(anmj), (!((size_t)(const void *)((anmi) + 1) - (size_t)(const
void *)(anmi) == 1) || __s1_len >= 4) && (!((size_t
)(const void *)((anmj) + 1) - (size_t)(const void *)(anmj) ==
1) || __s2_len >= 4)) ? __builtin_strcmp (anmi, anmj) : (
__builtin_constant_p (anmi) && ((size_t)(const void *
)((anmi) + 1) - (size_t)(const void *)(anmi) == 1) &&
(__s1_len = strlen (anmi), __s1_len < 4) ? (__builtin_constant_p
(anmj) && ((size_t)(const void *)((anmj) + 1) - (size_t
)(const void *)(anmj) == 1) ? __builtin_strcmp (anmi, anmj) :
(__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (anmj); int __result = (((const unsigned
char *) (const char *) (anmi))[0] - __s2[0]); if (__s1_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) (anmi))[1] - __s2[1]); if (__s1_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) (anmi))[2] - __s2[2]); if (__s1_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (anmi))[3] - __s2[3]); } } __result; }))) : (
__builtin_constant_p (anmj) && ((size_t)(const void *
)((anmj) + 1) - (size_t)(const void *)(anmj) == 1) &&
(__s2_len = strlen (anmj), __s2_len < 4) ? (__builtin_constant_p
(anmi) && ((size_t)(const void *)((anmi) + 1) - (size_t
)(const void *)(anmi) == 1) ? __builtin_strcmp (anmi, anmj) :
(- (__extension__ ({ const unsigned char *__s2 = (const unsigned
char *) (const char *) (anmi); int __result = (((const unsigned
char *) (const char *) (anmj))[0] - __s2[0]); if (__s2_len >
0 && __result == 0) { __result = (((const unsigned char
*) (const char *) (anmj))[1] - __s2[1]); if (__s2_len > 1
&& __result == 0) { __result = (((const unsigned char
*) (const char *) (anmj))[2] - __s2[2]); if (__s2_len > 2
&& __result == 0) __result = (((const unsigned char *
) (const char *) (anmj))[3] - __s2[3]); } } __result; })))) :
__builtin_strcmp (anmi, anmj)))); }) : strncmp (anmi, anmj, anmil
-1))) == 0 );
355 }
356 }
357 /* set index for this atom */
358 noe_index[i] = groupnr;
359 if (bMatch)
360 {
361 /* set index for next two matching atoms */
362 for (j = 1; j < 3; j++)
363 {
364 noe_index[i+j] = groupnr;
365 }
366 /* skip matching atoms */
367 i += 2;
368 }
369 }
370 groupnr++;
371 }
372 }
373 else
374 {
375 /* make index without looking for equivalent atoms */
376 for (i = 0; i < isize; i++)
377 {
378 noe_index[i] = i;
379 }
380 groupnr = isize;
381 }
382 noe_index[isize] = groupnr;
383
384 if (debug)
385 {
386 /* dump new names */
387 for (i = 0; i < isize; i++)
388 {
389 rnri = atoms->atom[index[i]].resind;
390 fprintf(debug, "%s %s %d -> %s\n", *atoms->atomname[index[i]],
391 *atoms->resinfo[rnri].name, rnri, nnm[i] ? nnm[i] : "");
392 }
393 }
394
395 for (i = 0; i < isize; i++)
396 {
397 gi = noe_index[i];
398 if (!noe_gr[gi].aname)
399 {
400 noe_gr[gi].ianr = i;
401 noe_gr[gi].anr = index[i];
402 if (nnm[i])
403 {
404 noe_gr[gi].aname = strdup(nnm[i])(__extension__ (__builtin_constant_p (nnm[i]) && ((size_t
)(const void *)((nnm[i]) + 1) - (size_t)(const void *)(nnm[i]
) == 1) ? (((const char *) (nnm[i]))[0] == '\0' ? (char *) calloc
((size_t) 1, (size_t) 1) : ({ size_t __len = strlen (nnm[i])
+ 1; char *__retval = (char *) malloc (__len); if (__retval !=
((void*)0)) __retval = (char *) memcpy (__retval, nnm[i], __len
); __retval; })) : __strdup (nnm[i])));
405 }
406 else
407 {
408 noe_gr[gi].aname = strdup(*atoms->atomname[index[i]])(__extension__ (__builtin_constant_p (*atoms->atomname[index
[i]]) && ((size_t)(const void *)((*atoms->atomname
[index[i]]) + 1) - (size_t)(const void *)(*atoms->atomname
[index[i]]) == 1) ? (((const char *) (*atoms->atomname[index
[i]]))[0] == '\0' ? (char *) calloc ((size_t) 1, (size_t) 1) :
({ size_t __len = strlen (*atoms->atomname[index[i]]) + 1
; char *__retval = (char *) malloc (__len); if (__retval != (
(void*)0)) __retval = (char *) memcpy (__retval, *atoms->atomname
[index[i]], __len); __retval; })) : __strdup (*atoms->atomname
[index[i]])));
409 if (noe_index[i] == noe_index[i+1])
410 {
411 noe_gr[gi].aname[strlen(noe_gr[gi].aname)-1] = '*';
412 }
413 }
414 noe_gr[gi].rnr = atoms->atom[index[i]].resind;
415 noe_gr[gi].rname = strdup(*atoms->resinfo[noe_gr[gi].rnr].name)(__extension__ (__builtin_constant_p (*atoms->resinfo[noe_gr
[gi].rnr].name) && ((size_t)(const void *)((*atoms->
resinfo[noe_gr[gi].rnr].name) + 1) - (size_t)(const void *)(*
atoms->resinfo[noe_gr[gi].rnr].name) == 1) ? (((const char
*) (*atoms->resinfo[noe_gr[gi].rnr].name))[0] == '\0' ? (
char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len = strlen
(*atoms->resinfo[noe_gr[gi].rnr].name) + 1; char *__retval
= (char *) malloc (__len); if (__retval != ((void*)0)) __retval
= (char *) memcpy (__retval, *atoms->resinfo[noe_gr[gi].rnr
].name, __len); __retval; })) : __strdup (*atoms->resinfo[
noe_gr[gi].rnr].name)));
416 /* dump group definitions */
417 if (debug)
418 {
419 fprintf(debug, "%d %d %d %d %s %s %d\n", i, gi,
420 noe_gr[gi].ianr, noe_gr[gi].anr, noe_gr[gi].aname,
421 noe_gr[gi].rname, noe_gr[gi].rnr);
422 }
423 }
424 }
425 for (i = 0; i < isize; i++)
426 {
427 sfree(nnm[i])save_free("nnm[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 427, (nnm[i]));
428 }
429 sfree(nnm)save_free("nnm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 429, (nnm));
430
431 return groupnr;
432}
433
434/* #define NSCALE 3 */
435/* static char *noe_scale[NSCALE+1] = { */
436/* "strong", "medium", "weak", "none" */
437/* }; */
438#define NSCALE6 6
439
440static char *noe2scale(real r3, real r6, real rmax)
441{
442 int i, s3, s6;
443 static char buf[NSCALE6+1];
444
445 /* r goes from 0 to rmax
446 NSCALE*r/rmax goes from 0 to NSCALE
447 NSCALE - NSCALE*r/rmax goes from NSCALE to 0 */
448 s3 = NSCALE6 - min(NSCALE, (int)(NSCALE*r3/rmax))(((6) < ((int)(6*r3/rmax))) ? (6) : ((int)(6*r3/rmax)) );
449 s6 = NSCALE6 - min(NSCALE, (int)(NSCALE*r6/rmax))(((6) < ((int)(6*r6/rmax))) ? (6) : ((int)(6*r6/rmax)) );
450
451 for (i = 0; i < s3; i++)
452 {
453 buf[i] = '=';
454 }
455 for (; i < s6; i++)
456 {
457 buf[i] = '-';
458 }
459 buf[i] = '\0';
460
461 return buf;
462}
463
464static void calc_noe(int isize, atom_id *noe_index,
465 real **dtot1_3, real **dtot1_6, int gnr, t_noe **noe)
466{
467 int i, j, gi, gj;
468
469 /* make half matrix of noe-group distances from atom distances */
470 for (i = 0; i < isize; i++)
471 {
472 gi = noe_index[i];
473 for (j = i; j < isize; j++)
474 {
475 gj = noe_index[j];
476 noe[gi][gj].nr++;
477 noe[gi][gj].i_3 += pow(dtot1_3[i][j], -3);
478 noe[gi][gj].i_6 += pow(dtot1_6[i][j], -6);
479 }
480 }
481
482 /* make averages */
483 for (i = 0; i < gnr; i++)
484 {
485 for (j = i+1; j < gnr; j++)
486 {
487 noe[i][j].r_3 = pow(noe[i][j].i_3/noe[i][j].nr, -1.0/3.0);
488 noe[i][j].r_6 = pow(noe[i][j].i_6/noe[i][j].nr, -1.0/6.0);
489 noe[j][i] = noe[i][j];
490 }
491 }
492}
493
494static void write_noe(FILE *fp, int gnr, t_noe **noe, t_noe_gr *noe_gr, real rmax)
495{
496 int i, j;
497 real r3, r6, min3, min6;
498 char buf[10], b3[10], b6[10];
499 t_noe_gr gri, grj;
500
501 min3 = min6 = 1e6;
502 fprintf(fp,
503 ";%4s %3s %4s %4s%3s %4s %4s %4s %4s%3s %5s %5s %8s %2s %2s %s\n",
504 "ianr", "anr", "anm", "rnm", "rnr", "ianr", "anr", "anm", "rnm", "rnr",
505 "1/r^3", "1/r^6", "intnsty", "Dr", "Da", "scale");
506 for (i = 0; i < gnr; i++)
507 {
508 gri = noe_gr[i];
509 for (j = i+1; j < gnr; j++)
510 {
511 grj = noe_gr[j];
512 r3 = noe[i][j].r_3;
513 r6 = noe[i][j].r_6;
514 min3 = min(r3, min3)(((r3) < (min3)) ? (r3) : (min3) );
515 min6 = min(r6, min6)(((r6) < (min6)) ? (r6) : (min6) );
516 if (r3 < rmax || r6 < rmax)
517 {
518 if (grj.rnr == gri.rnr)
519 {
520 sprintf(buf, "%2d", grj.anr-gri.anr);
521 }
522 else
523 {
524 buf[0] = '\0';
525 }
526 if (r3 < rmax)
527 {
528 sprintf(b3, "%-5.3f", r3);
529 }
530 else
531 {
532 strcpy(b3, "-");
533 }
534 if (r6 < rmax)
535 {
536 sprintf(b6, "%-5.3f", r6);
537 }
538 else
539 {
540 strcpy(b6, "-");
541 }
542 fprintf(fp,
543 "%4d %4d %4s %4s%3d %4d %4d %4s %4s%3d %5s %5s %8d %2d %2s %s\n",
544 gri.ianr+1, gri.anr+1, gri.aname, gri.rname, gri.rnr+1,
545 grj.ianr+1, grj.anr+1, grj.aname, grj.rname, grj.rnr+1,
546 b3, b6, (int)(noe[i][j].i_6+0.5), grj.rnr-gri.rnr, buf,
547 noe2scale(r3, r6, rmax));
548 }
549 }
550 }
551#define MINI ((i == 3) ? min3 : min6)
552 for (i = 3; i <= 6; i += 3)
553 {
554 if (MINI > rmax)
555 {
556 fprintf(stdoutstdout, "NOTE: no 1/r^%d averaged distances found below %g, "
557 "smallest was %g\n", i, rmax, MINI );
558 }
559 else
560 {
561 fprintf(stdoutstdout, "Smallest 1/r^%d averaged distance was %g\n", i, MINI );
562 }
563 }
564#undef MINI
565}
566
567static void calc_rms(int nind, int nframes,
568 real **dtot, real **dtot2,
569 real **rmsmat, real *rmsmax,
570 real **rmscmat, real *rmscmax,
571 real **meanmat, real *meanmax)
572{
573 int i, j;
574 real mean, mean2, rms, rmsc;
575/* N.B. dtot and dtot2 contain the total distance and the total squared
576 * distance respectively, BUT they return RMS and the scaled RMS resp.
577 */
578
579 *rmsmax = -1000;
580 *rmscmax = -1000;
581 *meanmax = -1000;
582
583 for (i = 0; (i < nind-1); i++)
584 {
585 for (j = i+1; (j < nind); j++)
586 {
587 mean = dtot[i][j]/nframes;
588 mean2 = dtot2[i][j]/nframes;
589 rms = sqrt(max(0, mean2-mean*mean)(((0) > (mean2-mean*mean)) ? (0) : (mean2-mean*mean) ));
590 rmsc = rms/mean;
591 if (mean > *meanmax)
592 {
593 *meanmax = mean;
594 }
595 if (rms > *rmsmax)
596 {
597 *rmsmax = rms;
598 }
599 if (rmsc > *rmscmax)
600 {
601 *rmscmax = rmsc;
602 }
603 meanmat[i][j] = meanmat[j][i] = mean;
604 rmsmat[i][j] = rmsmat[j][i] = rms;
605 rmscmat[i][j] = rmscmat[j][i] = rmsc;
606 }
607 }
608}
609
610real rms_diff(int natom, real **d, real **d_r)
611{
612 int i, j;
613 real r, r2;
614
615 r2 = 0.0;
616 for (i = 0; (i < natom-1); i++)
617 {
618 for (j = i+1; (j < natom); j++)
619 {
620 r = d[i][j]-d_r[i][j];
621 r2 += r*r;
622 }
623 }
624 r2 /= (natom*(natom-1))/2;
625
626 return sqrt(r2);
627}
628
629int gmx_rmsdist(int argc, char *argv[])
630{
631 const char *desc[] = {
632 "[THISMODULE] computes the root mean square deviation of atom distances,",
633 "which has the advantage that no fit is needed like in standard RMS",
634 "deviation as computed by [gmx-rms].",
635 "The reference structure is taken from the structure file.",
636 "The RMSD at time t is calculated as the RMS",
637 "of the differences in distance between atom-pairs in the reference",
638 "structure and the structure at time t.[PAR]",
639 "[THISMODULE] can also produce matrices of the rms distances, rms distances",
640 "scaled with the mean distance and the mean distances and matrices with",
641 "NMR averaged distances (1/r^3 and 1/r^6 averaging). Finally, lists",
642 "of atom pairs with 1/r^3 and 1/r^6 averaged distance below the",
643 "maximum distance ([TT]-max[tt], which will default to 0.6 in this case)",
644 "can be generated, by default averaging over equivalent hydrogens",
645 "(all triplets of hydrogens named *[123]). Additionally a list of",
646 "equivalent atoms can be supplied ([TT]-equiv[tt]), each line containing",
647 "a set of equivalent atoms specified as residue number and name and",
648 "atom name; e.g.:[PAR]",
649 "[TT]3 SER HB1 3 SER HB2[tt][PAR]",
650 "Residue and atom names must exactly match those in the structure",
651 "file, including case. Specifying non-sequential atoms is undefined."
652
653 };
654
655 int natom, i, j, teller, gi, gj;
656 real t;
657
658 t_topology top;
659 int ePBC;
660 t_atoms *atoms;
661 matrix box;
662 rvec *x;
663 FILE *fp;
664
665 t_trxstatus *status;
666 int isize, gnr = 0;
667 atom_id *index, *noe_index;
668 char *grpname;
669 real **d_r, **d, **dtot, **dtot2, **mean, **rms, **rmsc, *resnr;
670 real **dtot1_3 = NULL((void*)0), **dtot1_6 = NULL((void*)0);
671 real rmsnow, meanmax, rmsmax, rmscmax;
672 real max1_3, max1_6;
673 t_noe_gr *noe_gr = NULL((void*)0);
674 t_noe **noe = NULL((void*)0);
675 t_rgb rlo, rhi;
676 char buf[255];
677 gmx_bool bRMS, bScale, bMean, bNOE, bNMR3, bNMR6, bNMR;
678
679 static int nlevels = 40;
680 static real scalemax = -1.0;
681 static gmx_bool bSumH = TRUE1;
682 static gmx_bool bPBC = TRUE1;
683 output_env_t oenv;
684
685 t_pargs pa[] = {
686 { "-nlevels", FALSE0, etINT, {&nlevels},
687 "Discretize RMS in this number of levels" },
688 { "-max", FALSE0, etREAL, {&scalemax},
689 "Maximum level in matrices" },
690 { "-sumh", FALSE0, etBOOL, {&bSumH},
691 "Average distance over equivalent hydrogens" },
692 { "-pbc", FALSE0, etBOOL, {&bPBC},
693 "Use periodic boundary conditions when computing distances" }
694 };
695 t_filenm fnm[] = {
696 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
697 { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
698 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
699 { efDAT, "-equiv", "equiv", ffOPTRD(1<<1 | 1<<3) },
700 { efXVG, NULL((void*)0), "distrmsd", ffWRITE1<<2 },
701 { efXPM, "-rms", "rmsdist", ffOPTWR(1<<2| 1<<3) },
702 { efXPM, "-scl", "rmsscale", ffOPTWR(1<<2| 1<<3) },
703 { efXPM, "-mean", "rmsmean", ffOPTWR(1<<2| 1<<3) },
704 { efXPM, "-nmr3", "nmr3", ffOPTWR(1<<2| 1<<3) },
705 { efXPM, "-nmr6", "nmr6", ffOPTWR(1<<2| 1<<3) },
706 { efDAT, "-noe", "noe", ffOPTWR(1<<2| 1<<3) },
707 };
708#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
709
710 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW(1<<5) | PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_BE_NICE(1<<13),
711 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))
712 {
713 return 0;
714 }
715
716 bRMS = opt2bSet("-rms", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
717 bScale = opt2bSet("-scl", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
718 bMean = opt2bSet("-mean", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
719 bNOE = opt2bSet("-noe", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
720 bNMR3 = opt2bSet("-nmr3", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
721 bNMR6 = opt2bSet("-nmr6", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
722 bNMR = bNMR3 || bNMR6 || bNOE;
723
724 max1_3 = 0;
725 max1_6 = 0;
726
727 /* check input */
728 if (bNOE && scalemax < 0)
729 {
730 scalemax = 0.6;
731 fprintf(stderrstderr, "WARNING: using -noe without -max "
732 "makes no sense, setting -max to %g\n\n", scalemax);
733 }
734
735 /* get topology and index */
736 read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), buf, &top, &ePBC, &x, NULL((void*)0), box, FALSE0);
737
738 if (!bPBC)
739 {
740 ePBC = epbcNONE;
741 }
742 atoms = &(top.atoms);
743
744 get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &isize, &index, &grpname);
745
746 /* initialize arrays */
747 snew(d, isize)(d) = save_calloc("d", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 747, (isize), sizeof(*(d)));
748 snew(dtot, isize)(dtot) = save_calloc("dtot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 748, (isize), sizeof(*(dtot)));
749 snew(dtot2, isize)(dtot2) = save_calloc("dtot2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 749, (isize), sizeof(*(dtot2)));
750 if (bNMR)
751 {
752 snew(dtot1_3, isize)(dtot1_3) = save_calloc("dtot1_3", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 752, (isize), sizeof(*(dtot1_3)));
753 snew(dtot1_6, isize)(dtot1_6) = save_calloc("dtot1_6", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 753, (isize), sizeof(*(dtot1_6)));
754 }
755 snew(mean, isize)(mean) = save_calloc("mean", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 755, (isize), sizeof(*(mean)));
756 snew(rms, isize)(rms) = save_calloc("rms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 756, (isize), sizeof(*(rms)));
757 snew(rmsc, isize)(rmsc) = save_calloc("rmsc", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 757, (isize), sizeof(*(rmsc)));
758 snew(d_r, isize)(d_r) = save_calloc("d_r", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 758, (isize), sizeof(*(d_r)));
759 snew(resnr, isize)(resnr) = save_calloc("resnr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 759, (isize), sizeof(*(resnr)));
760 for (i = 0; (i < isize); i++)
761 {
762 snew(d[i], isize)(d[i]) = save_calloc("d[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 762, (isize), sizeof(*(d[i])));
763 snew(dtot[i], isize)(dtot[i]) = save_calloc("dtot[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 763, (isize), sizeof(*(dtot[i])));
764 snew(dtot2[i], isize)(dtot2[i]) = save_calloc("dtot2[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 764, (isize), sizeof(*(dtot2[i])));
765 if (bNMR)
766 {
767 snew(dtot1_3[i], isize)(dtot1_3[i]) = save_calloc("dtot1_3[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 767, (isize), sizeof(*(dtot1_3[i])));
768 snew(dtot1_6[i], isize)(dtot1_6[i]) = save_calloc("dtot1_6[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 768, (isize), sizeof(*(dtot1_6[i])));
769 }
770 snew(mean[i], isize)(mean[i]) = save_calloc("mean[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 770, (isize), sizeof(*(mean[i])));
771 snew(rms[i], isize)(rms[i]) = save_calloc("rms[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 771, (isize), sizeof(*(rms[i])));
772 snew(rmsc[i], isize)(rmsc[i]) = save_calloc("rmsc[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 772, (isize), sizeof(*(rmsc[i])));
773 snew(d_r[i], isize)(d_r[i]) = save_calloc("d_r[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 773, (isize), sizeof(*(d_r[i])));
774 resnr[i] = i+1;
775 }
776
777 /*set box type*/
778 calc_dist(isize, index, x, ePBC, box, d_r);
779 sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 779, (x));
780
781 /*open output files*/
782 fp = xvgropen(ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "RMS Deviation", "Time (ps)", "RMSD (nm)",
783 oenv);
784 if (output_env_get_print_xvgr_codes(oenv))
785 {
786 fprintf(fp, "@ subtitle \"of distances between %s atoms\"\n", grpname);
787 }
788
789 /*do a first step*/
790 natom = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box);
Value stored to 'natom' is never read
791
792 do
793 {
794 calc_dist_tot(isize, index, x, ePBC, box, d, dtot, dtot2, bNMR, dtot1_3, dtot1_6);
795
796 rmsnow = rms_diff(isize, d, d_r);
797 fprintf(fp, "%g %g\n", t, rmsnow);
798 }
799 while (read_next_x(oenv, status, &t, x, box));
800 fprintf(stderrstderr, "\n");
801
802 gmx_ffclose(fp);
803
804 teller = nframes_read(status);
805
806 close_trj(status);
807
808 calc_rms(isize, teller, dtot, dtot2, mean, &meanmax, rms, &rmsmax, rmsc, &rmscmax);
809 fprintf(stderrstderr, "rmsmax = %g, rmscmax = %g\n", rmsmax, rmscmax);
810
811 if (bNMR)
812 {
813 calc_nmr(isize, teller, dtot1_3, dtot1_6, &max1_3, &max1_6);
814 }
815
816 if (scalemax > -1.0)
817 {
818 rmsmax = scalemax;
819 rmscmax = scalemax;
820 meanmax = scalemax;
821 max1_3 = scalemax;
822 max1_6 = scalemax;
823 }
824
825 if (bNOE)
826 {
827 /* make list of noe atom groups */
828 snew(noe_index, isize+1)(noe_index) = save_calloc("noe_index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 828, (isize+1), sizeof(*(noe_index)));
829 snew(noe_gr, isize)(noe_gr) = save_calloc("noe_gr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 829, (isize), sizeof(*(noe_gr)));
830 gnr = analyze_noe_equivalent(opt2fn_null("-equiv", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
831 atoms, isize, index, bSumH, noe_index, noe_gr);
832 fprintf(stdoutstdout, "Found %d non-equivalent atom-groups in %d atoms\n",
833 gnr, isize);
834 /* make half matrix of of noe-group distances from atom distances */
835 snew(noe, gnr)(noe) = save_calloc("noe", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 835, (gnr), sizeof(*(noe)));
836 for (i = 0; i < gnr; i++)
837 {
838 snew(noe[i], gnr)(noe[i]) = save_calloc("noe[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rmsdist.c"
, 838, (gnr), sizeof(*(noe[i])));
839 }
840 calc_noe(isize, noe_index, dtot1_3, dtot1_6, gnr, noe);
841 }
842
843 rlo.r = 1.0, rlo.g = 1.0, rlo.b = 1.0;
844 rhi.r = 0.0, rhi.g = 0.0, rhi.b = 0.0;
845
846 if (bRMS)
847 {
848 write_xpm(opt2FILE("-rms", NFILE, fnm, "w")gmx_ffopen(opt2fn("-rms", ((int)(sizeof(fnm)/sizeof((fnm)[0])
)), fnm), "w"), 0,
849 "RMS of distance", "RMS (nm)", "Atom Index", "Atom Index",
850 isize, isize, resnr, resnr, rms, 0.0, rmsmax, rlo, rhi, &nlevels);
851 }
852
853 if (bScale)
854 {
855 write_xpm(opt2FILE("-scl", NFILE, fnm, "w")gmx_ffopen(opt2fn("-scl", ((int)(sizeof(fnm)/sizeof((fnm)[0])
)), fnm), "w"), 0,
856 "Relative RMS", "RMS", "Atom Index", "Atom Index",
857 isize, isize, resnr, resnr, rmsc, 0.0, rmscmax, rlo, rhi, &nlevels);
858 }
859
860 if (bMean)
861 {
862 write_xpm(opt2FILE("-mean", NFILE, fnm, "w")gmx_ffopen(opt2fn("-mean", ((int)(sizeof(fnm)/sizeof((fnm)[0]
))), fnm), "w"), 0,
863 "Mean Distance", "Distance (nm)", "Atom Index", "Atom Index",
864 isize, isize, resnr, resnr, mean, 0.0, meanmax, rlo, rhi, &nlevels);
865 }
866
867 if (bNMR3)
868 {
869 write_xpm(opt2FILE("-nmr3", NFILE, fnm, "w")gmx_ffopen(opt2fn("-nmr3", ((int)(sizeof(fnm)/sizeof((fnm)[0]
))), fnm), "w"), 0, "1/r^3 averaged distances",
870 "Distance (nm)", "Atom Index", "Atom Index",
871 isize, isize, resnr, resnr, dtot1_3, 0.0, max1_3, rlo, rhi, &nlevels);
872 }
873 if (bNMR6)
874 {
875 write_xpm(opt2FILE("-nmr6", NFILE, fnm, "w")gmx_ffopen(opt2fn("-nmr6", ((int)(sizeof(fnm)/sizeof((fnm)[0]
))), fnm), "w"), 0, "1/r^6 averaged distances",
876 "Distance (nm)", "Atom Index", "Atom Index",
877 isize, isize, resnr, resnr, dtot1_6, 0.0, max1_6, rlo, rhi, &nlevels);
878 }
879
880 if (bNOE)
881 {
882 write_noe(opt2FILE("-noe", NFILE, fnm, "w")gmx_ffopen(opt2fn("-noe", ((int)(sizeof(fnm)/sizeof((fnm)[0])
)), fnm), "w"), gnr, noe, noe_gr, scalemax);
883 }
884
885 do_view(oenv, ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
886
887 return 0;
888}