Bug Summary

File:gromacs/gmxana/gmx_anadock.c
Location:line 246, column 9
Description:Value stored to 'dist' 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) 2012,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 <stdlib.h>
42#include <string.h>
43
44#include "gromacs/fileio/confio.h"
45#include "copyrite.h"
46#include "gromacs/utility/futil.h"
47#include "gromacs/utility/fatalerror.h"
48#include "gromacs/utility/smalloc.h"
49#include "gromacs/utility/cstringutil.h"
50#include "gromacs/math/vec.h"
51#include "gromacs/statistics/statistics.h"
52#include "gromacs/commandline/pargs.h"
53#include "typedefs.h"
54#include "gromacs/fileio/xvgr.h"
55#include "macros.h"
56
57static const char *etitles[] = { "E-docked", "Free Energy" };
58
59typedef struct {
60 real edocked, efree;
61 int index, cluster_id;
62 t_atoms atoms;
63 rvec *x;
64 int ePBC;
65 matrix box;
66} t_pdbfile;
67
68static t_pdbfile *read_pdbf(const char *fn)
69{
70 t_pdbfile *pdbf;
71 double e;
72 char buf[256], *ptr;
73 int natoms;
74 FILE *fp;
75
76 snew(pdbf, 1)(pdbf) = save_calloc("pdbf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 76, (1), sizeof(*(pdbf)));
77 get_stx_coordnum (fn, &natoms);
78 init_t_atoms(&(pdbf->atoms), natoms, FALSE0);
79 snew(pdbf->x, natoms)(pdbf->x) = save_calloc("pdbf->x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 79, (natoms), sizeof(*(pdbf->x)));
80 read_stx_conf(fn, buf, &pdbf->atoms, pdbf->x, NULL((void*)0), &pdbf->ePBC, pdbf->box);
81 fp = gmx_ffopen(fn, "r");
82 do
83 {
84 ptr = fgets2(buf, 255, fp);
85 if (ptr)
86 {
87 if (strstr(buf, "Intermolecular") != NULL((void*)0))
88 {
89 ptr = strchr(buf, '=')(__extension__ (__builtin_constant_p ('=') && !__builtin_constant_p
(buf) && ('=') == '\0' ? (char *) __rawmemchr (buf, '='
) : __builtin_strchr (buf, '=')));
90 sscanf(ptr+1, "%lf", &e);
91 pdbf->edocked = e;
92 }
93 else if (strstr(buf, "Estimated Free") != NULL((void*)0))
94 {
95 ptr = strchr(buf, '=')(__extension__ (__builtin_constant_p ('=') && !__builtin_constant_p
(buf) && ('=') == '\0' ? (char *) __rawmemchr (buf, '='
) : __builtin_strchr (buf, '=')));
96 sscanf(ptr+1, "%lf", &e);
97 pdbf->efree = e;
98 }
99 }
100 }
101 while (ptr != NULL((void*)0));
102 gmx_ffclose(fp);
103
104 return pdbf;
105}
106
107static t_pdbfile **read_em_all(const char *fn, int *npdbf)
108{
109 t_pdbfile **pdbf = 0;
110 int i, maxpdbf;
111 char buf[256], name[256];
112 gmx_bool bExist;
113
114 strcpy(buf, fn);
115 buf[strlen(buf)-4] = '\0';
116 maxpdbf = 100;
117 snew(pdbf, maxpdbf)(pdbf) = save_calloc("pdbf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 117, (maxpdbf), sizeof(*(pdbf)));
118 i = 0;
119 do
120 {
121 sprintf(name, "%s_%d.pdb", buf, i+1);
122 if ((bExist = gmx_fexist(name)) == TRUE1)
123 {
124 pdbf[i] = read_pdbf(name);
125 pdbf[i]->index = i+1;
126 i++;
127 if (i >= maxpdbf)
128 {
129 maxpdbf += 100;
130 srenew(pdbf, maxpdbf)(pdbf) = save_realloc("pdbf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 130, (pdbf), (maxpdbf), sizeof(*(pdbf)));
131 }
132 }
133 }
134 while (bExist);
135
136 *npdbf = i;
137
138 printf("Found %d pdb files\n", i);
139
140 return pdbf;
141}
142
143static gmx_bool bFreeSort = FALSE0;
144
145static int pdbf_comp(const void *a, const void *b)
146{
147 t_pdbfile *pa, *pb;
148 real x;
149 int dc;
150
151 pa = *(t_pdbfile **)a;
152 pb = *(t_pdbfile **)b;
153
154 dc = pa->cluster_id - pb->cluster_id;
155
156 if (dc == 0)
157 {
158 if (bFreeSort)
159 {
160 x = pa->efree - pb->efree;
161 }
162 else
163 {
164 x = pa->edocked - pb->edocked;
165 }
166
167 if (x < 0)
168 {
169 return -1;
170 }
171 else if (x > 0)
172 {
173 return 1;
174 }
175 else
176 {
177 return 0;
178 }
179 }
180 else
181 {
182 return dc;
183 }
184}
185
186static void analyse_em_all(int npdb, t_pdbfile *pdbf[], const char *edocked,
187 const char *efree, const output_env_t oenv)
188{
189 FILE *fp;
190 int i;
191
192 for (bFreeSort = FALSE0; (bFreeSort <= TRUE1); bFreeSort++)
193 {
194 qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
195 fp = xvgropen(bFreeSort ? efree : edocked,
196 etitles[bFreeSort], "()", "E (kJ/mol)", oenv);
197 for (i = 0; (i < npdb); i++)
198 {
199 fprintf(fp, "%12lf\n", bFreeSort ? pdbf[i]->efree : pdbf[i]->edocked);
200 }
201 gmx_ffclose(fp);
202 }
203}
204
205static void clust_stat(FILE *fp, int start, int end, t_pdbfile *pdbf[])
206{
207 int i;
208 gmx_stats_t ed, ef;
209 real aver, sigma;
210
211 ed = gmx_stats_init();
212 ef = gmx_stats_init();
213 for (i = start; (i < end); i++)
214 {
215 gmx_stats_add_point(ed, i-start, pdbf[i]->edocked, 0, 0);
216 gmx_stats_add_point(ef, i-start, pdbf[i]->efree, 0, 0);
217 }
218 gmx_stats_get_ase(ed, &aver, &sigma, NULL((void*)0));
219 fprintf(fp, " <%12s> = %8.3f (+/- %6.3f)\n", etitles[FALSE0], aver, sigma);
220 gmx_stats_get_ase(ef, &aver, &sigma, NULL((void*)0));
221 fprintf(fp, " <%12s> = %8.3f (+/- %6.3f)\n", etitles[TRUE1], aver, sigma);
222 gmx_stats_done(ed);
223 gmx_stats_done(ef);
224 sfree(ed)save_free("ed", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 224, (ed));
225 sfree(ef)save_free("ef", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 225, (ef));
226}
227
228static real rmsd_dist(t_pdbfile *pa, t_pdbfile *pb, gmx_bool bRMSD)
229{
230 int i;
231 real rmsd, dist;
232 rvec acm, bcm, dx;
233
234 if (bRMSD)
235 {
236 rmsd = 0;
237 for (i = 0; (i < pa->atoms.nr); i++)
238 {
239 rvec_sub(pa->x[i], pb->x[i], dx);
240 rmsd += iprod(dx, dx);
241 }
242 rmsd = sqrt(rmsd/pa->atoms.nr);
243 }
244 else
245 {
246 dist = 0;
Value stored to 'dist' is never read
247 clear_rvec(acm);
248 clear_rvec(bcm);
249 for (i = 0; (i < pa->atoms.nr); i++)
250 {
251 rvec_inc(acm, pa->x[i]);
252 rvec_inc(bcm, pb->x[i]);
253 }
254 rvec_sub(acm, bcm, dx);
255 for (i = 0; (i < DIM3); i++)
256 {
257 dx[i] /= pa->atoms.nr;
258 }
259 rmsd = norm(dx);
260 }
261 return rmsd;
262}
263
264static void line(FILE *fp)
265{
266 fprintf(fp, " - - - - - - - - - -\n");
267}
268
269static void cluster_em_all(FILE *fp, int npdb, t_pdbfile *pdbf[],
270 gmx_bool bFree, gmx_bool bRMSD, real cutoff)
271{
272 int i, j, k;
273 int *cndx, ncluster;
274 real rmsd;
275
276 bFreeSort = bFree;
277 qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
278
279 fprintf(fp, "Statistics over all structures:\n");
280 clust_stat(fp, 0, npdb, pdbf);
281 line(fp);
282
283 /* Index to first structure in a cluster */
284 snew(cndx, npdb)(cndx) = save_calloc("cndx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 284, (npdb), sizeof(*(cndx)));
285 ncluster = 1;
286
287 for (i = 1; (i < npdb); i++)
288 {
289 for (j = 0; (j < ncluster); j++)
290 {
291 rmsd = rmsd_dist(pdbf[cndx[j]], pdbf[i], bRMSD);
292 if (rmsd <= cutoff)
293 {
294 /* Structure i is in cluster j */
295 pdbf[i]->cluster_id = pdbf[cndx[j]]->cluster_id;
296 break;
297 }
298 }
299 if (j == ncluster)
300 {
301 /* New cluster! Cool! */
302 cndx[ncluster] = i;
303 pdbf[i]->cluster_id = ncluster;
304 ncluster++;
305 }
306 }
307 cndx[ncluster] = npdb;
308 qsort(pdbf, npdb, sizeof(pdbf[0]), pdbf_comp);
309
310 j = 0;
311 cndx[j++] = 0;
312 for (i = 1; (i < npdb); i++)
313 {
314 if (pdbf[i]->cluster_id != pdbf[i-1]->cluster_id)
315 {
316 cndx[j++] = i;
317 }
318 }
319 cndx[j] = npdb;
320 if (j != ncluster)
321 {
322 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 322, "Consistency error: j = %d, ncluster = %d", j, ncluster);
323 }
324
325 fprintf(fp, "I found %d clusters based on %s and %s with a %.3f nm cut-off\n",
326 ncluster, etitles[bFree], bRMSD ? "RMSD" : "distance", cutoff);
327 line(fp);
328 for (j = 0; (j < ncluster); j++)
329 {
330 fprintf(fp, "Cluster: %3d %s: %10.5f kJ/mol %3d elements\n",
331 j, etitles[bFree],
332 bFree ? pdbf[cndx[j]]->efree : pdbf[cndx[j]]->edocked,
333 cndx[j+1]-cndx[j]);
334 clust_stat(fp, cndx[j], cndx[j+1], pdbf);
335 for (k = cndx[j]; (k < cndx[j+1]); k++)
336 {
337 fprintf(fp, " %3d", pdbf[k]->index);
338 }
339 fprintf(fp, "\n");
340 line(fp);
341 }
342 sfree(cndx)save_free("cndx", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_anadock.c"
, 342, (cndx));
343}
344
345int gmx_anadock(int argc, char *argv[])
346{
347 const char *desc[] = {
348 "[THISMODULE] analyses the results of an Autodock run and clusters the",
349 "structures together, based on distance or RMSD. The docked energy",
350 "and free energy estimates are analysed, and for each cluster the",
351 "energy statistics are printed.[PAR]",
352 "An alternative approach to this is to cluster the structures first",
353 "using [gmx-cluster] and then sort the clusters on either lowest",
354 "energy or average energy."
355 };
356 t_filenm fnm[] = {
357 { efPDB, "-f", NULL((void*)0), ffREAD1<<1 },
358 { efXVG, "-od", "edocked", ffWRITE1<<2 },
359 { efXVG, "-of", "efree", ffWRITE1<<2 },
360 { efLOG, "-g", "anadock", ffWRITE1<<2 }
361 };
362 output_env_t oenv;
363#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
364 static gmx_bool bFree = FALSE0, bRMS = TRUE1;
365 static real cutoff = 0.2;
366 t_pargs pa[] = {
367 { "-free", FALSE0, etBOOL, {&bFree},
368 "Use Free energy estimate from autodock for sorting the classes" },
369 { "-rms", FALSE0, etBOOL, {&bRMS},
370 "Cluster on RMS or distance" },
371 { "-cutoff", FALSE0, etREAL, {&cutoff},
372 "Maximum RMSD/distance for belonging to the same cluster" }
373 };
374#define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])))
375
376 FILE *fp;
377 t_pdbfile **pdbf = NULL((void*)0);
378 int npdbf;
379
380 if (!parse_common_args(&argc, argv, 0, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, NPA((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0,
381 NULL((void*)0), &oenv))
382 {
383 return 0;
384 }
385
386 fp = gmx_ffopen(opt2fn("-g", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "w");
387 please_cite(stdoutstdout, "Hetenyi2002b");
388 please_cite(fp, "Hetenyi2002b");
389
390 pdbf = read_em_all(opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &npdbf);
391
392 analyse_em_all(npdbf, pdbf, opt2fn("-od", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), opt2fn("-of", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm),
393 oenv);
394
395 cluster_em_all(fp, npdbf, pdbf, bFree, bRMS, cutoff);
396
397 gmx_thanx(fp);
398 gmx_ffclose(fp);
399
400 return 0;
401}