Bug Summary

File:gromacs/gmxana/gmx_trjorder.c
Location:line 326, column 21
Description:Value stored to 'sa' 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 <stdlib.h>
43#include <string.h>
44
45#include "gromacs/commandline/pargs.h"
46#include "typedefs.h"
47#include "gromacs/utility/smalloc.h"
48#include "macros.h"
49#include "gromacs/math/vec.h"
50#include "pbc.h"
51#include "gromacs/utility/futil.h"
52#include "index.h"
53#include "mshift.h"
54#include "gromacs/fileio/xvgr.h"
55#include "princ.h"
56#include "rmpbc.h"
57#include "txtdump.h"
58#include "gromacs/fileio/tpxio.h"
59#include "gromacs/fileio/trxio.h"
60#include "gmx_ana.h"
61
62#include "gromacs/utility/fatalerror.h"
63
64typedef struct {
65 atom_id i;
66 real d2;
67} t_order;
68
69t_order *order;
70
71static int ocomp(const void *a, const void *b)
72{
73 t_order *oa, *ob;
74
75 oa = (t_order *)a;
76 ob = (t_order *)b;
77
78 if (oa->d2 < ob->d2)
79 {
80 return -1;
81 }
82 else
83 {
84 return 1;
85 }
86}
87
88int gmx_trjorder(int argc, char *argv[])
89{
90 const char *desc[] = {
91 "[THISMODULE] orders molecules according to the smallest distance",
92 "to atoms in a reference group",
93 "or on z-coordinate (with option [TT]-z[tt]).",
94 "With distance ordering, it will ask for a group of reference",
95 "atoms and a group of molecules. For each frame of the trajectory",
96 "the selected molecules will be reordered according to the shortest",
97 "distance between atom number [TT]-da[tt] in the molecule and all the",
98 "atoms in the reference group. The center of mass of the molecules can",
99 "be used instead of a reference atom by setting [TT]-da[tt] to 0.",
100 "All atoms in the trajectory are written",
101 "to the output trajectory.[PAR]",
102 "[THISMODULE] can be useful for e.g. analyzing the n waters closest to a",
103 "protein.",
104 "In that case the reference group would be the protein and the group",
105 "of molecules would consist of all the water atoms. When an index group",
106 "of the first n waters is made, the ordered trajectory can be used",
107 "with any Gromacs program to analyze the n closest waters.",
108 "[PAR]",
109 "If the output file is a [TT].pdb[tt] file, the distance to the reference target",
110 "will be stored in the B-factor field in order to color with e.g. Rasmol.",
111 "[PAR]",
112 "With option [TT]-nshell[tt] the number of molecules within a shell",
113 "of radius [TT]-r[tt] around the reference group are printed."
114 };
115 static int na = 3, ref_a = 1;
116 static real rcut = 0;
117 static gmx_bool bCOM = FALSE0, bZ = FALSE0;
118 t_pargs pa[] = {
119 { "-na", FALSE0, etINT, {&na},
120 "Number of atoms in a molecule" },
121 { "-da", FALSE0, etINT, {&ref_a},
122 "Atom used for the distance calculation, 0 is COM" },
123 { "-com", FALSE0, etBOOL, {&bCOM},
124 "Use the distance to the center of mass of the reference group" },
125 { "-r", FALSE0, etREAL, {&rcut},
126 "Cutoff used for the distance calculation when computing the number of molecules in a shell around e.g. a protein" },
127 { "-z", FALSE0, etBOOL, {&bZ},
128 "Order molecules on z-coordinate" }
129 };
130 FILE *fp;
131 t_trxstatus *out;
132 t_trxstatus *status;
133 gmx_bool bNShell, bPDBout;
134 t_topology top;
135 int ePBC;
136 rvec *x, *xsol, xcom, dx;
137 matrix box;
138 t_pbc pbc;
139 gmx_rmpbc_t gpbc;
140 real t, totmass, mass, rcut2 = 0, n2;
141 int natoms, nwat, ncut;
142 char **grpname, title[256];
143 int i, j, d, *isize, isize_ref = 0, isize_sol;
144 atom_id sa, sr, *swi, **index, *ind_ref = NULL((void*)0), *ind_sol;
145 output_env_t oenv;
146 t_filenm fnm[] = {
147 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
148 { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
149 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
150 { efTRO, "-o", "ordered", ffOPTWR(1<<2| 1<<3) },
151 { efXVG, "-nshell", "nshell", ffOPTWR(1<<2| 1<<3) }
152 };
153#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
154
155 if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_BE_NICE(1<<13),
156 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))
157 {
158 return 0;
159 }
160
161 read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &x, NULL((void*)0), box, TRUE1);
162 sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 162, (x));
163
164 /* get index groups */
165 printf("Select %sa group of molecules to be ordered:\n",
166 bZ ? "" : "a group of reference atoms and ");
167 snew(grpname, 2)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 167, (2), sizeof(*(grpname)));
168 snew(index, 2)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 168, (2), sizeof(*(index)));
169 snew(isize, 2)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 169, (2), sizeof(*(isize)));
170 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), bZ ? 1 : 2,
171 isize, index, grpname);
172
173 if (!bZ)
174 {
175 isize_ref = isize[0];
176 isize_sol = isize[1];
177 ind_ref = index[0];
178 ind_sol = index[1];
179 }
180 else
181 {
182 isize_sol = isize[0];
183 ind_sol = index[0];
184 }
185
186 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box);
187 if (natoms > top.atoms.nr)
188 {
189 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 189, "Number of atoms in the run input file is larger than in the trjactory");
190 }
191 for (i = 0; (i < 2); i++)
192 {
193 for (j = 0; (j < isize[i]); j++)
194 {
195 if (index[i][j] > natoms)
196 {
197 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 197, "An atom number in group %s is larger than the number of atoms in the trajectory");
198 }
199 }
200 }
201
202 if ((isize_sol % na) != 0)
203 {
204 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 204, "Number of atoms in the molecule group (%d) is not a multiple of na (%d)",
205 isize[1], na);
206 }
207
208 nwat = isize_sol/na;
209 if (ref_a > na)
210 {
211 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 211, "The reference atom can not be larger than the number of atoms in a molecule");
212 }
213 ref_a--;
214 snew(xsol, nwat)(xsol) = save_calloc("xsol", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 214, (nwat), sizeof(*(xsol)));
215 snew(order, nwat)(order) = save_calloc("order", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 215, (nwat), sizeof(*(order)));
216 snew(swi, natoms)(swi) = save_calloc("swi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 216, (natoms), sizeof(*(swi)));
217 for (i = 0; (i < natoms); i++)
218 {
219 swi[i] = i;
220 }
221
222 out = NULL((void*)0);
223 fp = NULL((void*)0);
224 bNShell = ((opt2bSet("-nshell", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) ||
225 (opt2parg_bSet("-r", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa)));
226 bPDBout = FALSE0;
227 if (bNShell)
228 {
229 rcut2 = rcut*rcut;
230 fp = xvgropen(opt2fn("-nshell", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Number of molecules",
231 "Time (ps)", "N", oenv);
232 printf("Will compute the number of molecules within a radius of %g\n",
233 rcut);
234 }
235 if (!bNShell || opt2bSet("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm))
236 {
237 bPDBout = (fn2ftp(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) == efPDB);
238 if (bPDBout && !top.atoms.pdbinfo)
239 {
240 fprintf(stderrstderr, "Creating pdbfino records\n");
241 snew(top.atoms.pdbinfo, top.atoms.nr)(top.atoms.pdbinfo) = save_calloc("top.atoms.pdbinfo", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_trjorder.c"
, 241, (top.atoms.nr), sizeof(*(top.atoms.pdbinfo)));
242 }
243 out = open_trx(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "w");
244 }
245 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
246 do
247 {
248 gmx_rmpbc(gpbc, natoms, box, x);
249 set_pbc(&pbc, ePBC, box);
250
251 if (ref_a == -1)
252 {
253 /* Calculate the COM of all solvent molecules */
254 for (i = 0; i < nwat; i++)
255 {
256 totmass = 0;
257 clear_rvec(xsol[i]);
258 for (j = 0; j < na; j++)
259 {
260 sa = ind_sol[i*na+j];
261 mass = top.atoms.atom[sa].m;
262 totmass += mass;
263 for (d = 0; d < DIM3; d++)
264 {
265 xsol[i][d] += mass*x[sa][d];
266 }
267 }
268 svmul(1/totmass, xsol[i], xsol[i]);
269 }
270 }
271 else
272 {
273 /* Copy the reference atom of all solvent molecules */
274 for (i = 0; i < nwat; i++)
275 {
276 copy_rvec(x[ind_sol[i*na+ref_a]], xsol[i]);
277 }
278 }
279
280 if (bZ)
281 {
282 for (i = 0; (i < nwat); i++)
283 {
284 sa = ind_sol[na*i];
285 order[i].i = sa;
286 order[i].d2 = xsol[i][ZZ2];
287 }
288 }
289 else if (bCOM)
290 {
291 totmass = 0;
292 clear_rvec(xcom);
293 for (i = 0; i < isize_ref; i++)
294 {
295 mass = top.atoms.atom[ind_ref[i]].m;
296 totmass += mass;
297 for (j = 0; j < DIM3; j++)
298 {
299 xcom[j] += mass*x[ind_ref[i]][j];
300 }
301 }
302 svmul(1/totmass, xcom, xcom);
303 for (i = 0; (i < nwat); i++)
304 {
305 sa = ind_sol[na*i];
306 pbc_dx(&pbc, xcom, xsol[i], dx);
307 order[i].i = sa;
308 order[i].d2 = norm2(dx);
309 }
310 }
311 else
312 {
313 /* Set distance to first atom */
314 for (i = 0; (i < nwat); i++)
315 {
316 sa = ind_sol[na*i];
317 pbc_dx(&pbc, x[ind_ref[0]], xsol[i], dx);
318 order[i].i = sa;
319 order[i].d2 = norm2(dx);
320 }
321 for (j = 1; (j < isize_ref); j++)
322 {
323 sr = ind_ref[j];
324 for (i = 0; (i < nwat); i++)
325 {
326 sa = ind_sol[na*i];
Value stored to 'sa' is never read
327 pbc_dx(&pbc, x[sr], xsol[i], dx);
328 n2 = norm2(dx);
329 if (n2 < order[i].d2)
330 {
331 order[i].d2 = n2;
332 }
333 }
334 }
335 }
336
337 if (bNShell)
338 {
339 ncut = 0;
340 for (i = 0; (i < nwat); i++)
341 {
342 if (order[i].d2 <= rcut2)
343 {
344 ncut++;
345 }
346 }
347 fprintf(fp, "%10.3f %8d\n", t, ncut);
348 }
349 if (out)
350 {
351 qsort(order, nwat, sizeof(*order), ocomp);
352 for (i = 0; (i < nwat); i++)
353 {
354 for (j = 0; (j < na); j++)
355 {
356 swi[ind_sol[na*i]+j] = order[i].i+j;
357 }
358 }
359
360 /* Store the distance as the B-factor */
361 if (bPDBout)
362 {
363 for (i = 0; (i < nwat); i++)
364 {
365 for (j = 0; (j < na); j++)
366 {
367 top.atoms.pdbinfo[order[i].i+j].bfac = sqrt(order[i].d2);
368 }
369 }
370 }
371 write_trx(out, natoms, swi, &top.atoms, 0, t, box, x, NULL((void*)0), NULL((void*)0));
372 }
373 }
374 while (read_next_x(oenv, status, &t, x, box));
375 close_trj(status);
376 if (out)
377 {
378 close_trx(out);
379 }
380 if (fp)
381 {
382 gmx_ffclose(fp);
383 }
384 gmx_rmpbc_done(gpbc);
385
386 return 0;
387}