Bug Summary

File:gromacs/gmxana/gmx_sorient.c
Location:line 294, column 9
Description:Value stored to 'outp' 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 "macros.h"
42#include "gromacs/commandline/pargs.h"
43#include "gromacs/utility/smalloc.h"
44#include "gstat.h"
45#include "gromacs/math/vec.h"
46#include "gromacs/fileio/xvgr.h"
47#include "viewit.h"
48#include "pbc.h"
49#include "index.h"
50#include "gromacs/fileio/tpxio.h"
51#include "gromacs/fileio/trxio.h"
52#include "gmx_ana.h"
53
54#include "gromacs/utility/fatalerror.h"
55
56static void calc_com_pbc(int nrefat, t_topology *top, rvec x[], t_pbc *pbc,
57 atom_id index[], rvec xref, gmx_bool bPBC)
58{
59 const real tol = 1e-4;
60 gmx_bool bChanged;
61 int m, j, ai, iter;
62 real mass, mtot;
63 rvec dx, xtest;
64
65 /* First simple calculation */
66 clear_rvec(xref);
67 mtot = 0;
68 for (m = 0; (m < nrefat); m++)
69 {
70 ai = index[m];
71 mass = top->atoms.atom[ai].m;
72 for (j = 0; (j < DIM3); j++)
73 {
74 xref[j] += mass*x[ai][j];
75 }
76 mtot += mass;
77 }
78 svmul(1/mtot, xref, xref);
79 /* Now check if any atom is more than half the box from the COM */
80 if (bPBC)
81 {
82 iter = 0;
83 do
84 {
85 bChanged = FALSE0;
86 for (m = 0; (m < nrefat); m++)
87 {
88 ai = index[m];
89 mass = top->atoms.atom[ai].m/mtot;
90 pbc_dx(pbc, x[ai], xref, dx);
91 rvec_add(xref, dx, xtest);
92 for (j = 0; (j < DIM3); j++)
93 {
94 if (fabs(xtest[j]-x[ai][j]) > tol)
95 {
96 /* Here we have used the wrong image for contributing to the COM */
97 xref[j] += mass*(xtest[j]-x[ai][j]);
98 x[ai][j] = xtest[j];
99 bChanged = TRUE1;
100 }
101 }
102 }
103 if (bChanged)
104 {
105 printf("COM: %8.3f %8.3f %8.3f iter = %d\n", xref[XX0], xref[YY1], xref[ZZ2], iter);
106 }
107 iter++;
108 }
109 while (bChanged);
110 }
111}
112
113int gmx_sorient(int argc, char *argv[])
114{
115 t_topology top;
116 int ePBC = -1;
117 char title[STRLEN4096];
118 t_trxstatus *status;
119 int natoms;
120 real t;
121 rvec *xtop, *x;
122 matrix box;
123
124 FILE *fp;
125 int i, j, p, sa0, sa1, sa2, n, ntot, nf, m, *hist1, *hist2, *histn, nbin1, nbin2, nrbin;
126 real *histi1, *histi2, invbw, invrbw;
127 double sum1, sum2;
128 int *isize, nrefgrp, nrefat;
129 atom_id **index;
130 char **grpname;
131 real inp, outp, two_pi, nav, normfac, rmin2, rmax2, rcut, rcut2, r2, r, mass, mtot;
132 real c1, c2;
133 char str[STRLEN4096];
134 gmx_bool bTPS;
135 rvec xref, dx, dxh1, dxh2, outer;
136 gmx_rmpbc_t gpbc = NULL((void*)0);
137 t_pbc pbc;
138 const char *legr[] = {
139 "<cos(\\8q\\4\\s1\\N)>",
140 "<3cos\\S2\\N(\\8q\\4\\s2\\N)-1>"
141 };
142 const char *legc[] = {
143 "cos(\\8q\\4\\s1\\N)",
144 "3cos\\S2\\N(\\8q\\4\\s2\\N)-1"
145 };
146
147 const char *desc[] = {
148 "[THISMODULE] analyzes solvent orientation around solutes.",
149 "It calculates two angles between the vector from one or more",
150 "reference positions to the first atom of each solvent molecule:[PAR]",
151 "[GRK]theta[grk][SUB]1[sub]: the angle with the vector from the first atom of the solvent",
152 "molecule to the midpoint between atoms 2 and 3.[BR]",
153 "[GRK]theta[grk][SUB]2[sub]: the angle with the normal of the solvent plane, defined by the",
154 "same three atoms, or, when the option [TT]-v23[tt] is set, ",
155 "the angle with the vector between atoms 2 and 3.[PAR]",
156 "The reference can be a set of atoms or",
157 "the center of mass of a set of atoms. The group of solvent atoms should",
158 "consist of 3 atoms per solvent molecule.",
159 "Only solvent molecules between [TT]-rmin[tt] and [TT]-rmax[tt] are",
160 "considered for [TT]-o[tt] and [TT]-no[tt] each frame.[PAR]",
161 "[TT]-o[tt]: distribtion of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] for rmin<=r<=rmax.[PAR]",
162 "[TT]-no[tt]: distribution of [MATH][COS][GRK]theta[grk][SUB]2[sub][cos][math] for rmin<=r<=rmax.[PAR]",
163 "[TT]-ro[tt]: [MATH][CHEVRON][COS][GRK]theta[grk][SUB]1[sub][cos][chevron][math] and [MATH][CHEVRON]3[COS]^2[GRK]theta[grk][SUB]2[sub][cos]-1[chevron][math] as a function of the",
164 "distance.[PAR]",
165 "[TT]-co[tt]: the sum over all solvent molecules within distance r",
166 "of [MATH][COS][GRK]theta[grk][SUB]1[sub][cos][math] and [MATH]3[COS]^2([GRK]theta[grk][SUB]2[sub])-1[cos][math] as a function of r.[PAR]",
167 "[TT]-rc[tt]: the distribution of the solvent molecules as a function of r"
168 };
169
170 output_env_t oenv;
171 static gmx_bool bCom = FALSE0, bVec23 = FALSE0, bPBC = FALSE0;
172 static real rmin = 0.0, rmax = 0.5, binwidth = 0.02, rbinw = 0.02;
173 t_pargs pa[] = {
174 { "-com", FALSE0, etBOOL, {&bCom},
175 "Use the center of mass as the reference postion" },
176 { "-v23", FALSE0, etBOOL, {&bVec23},
177 "Use the vector between atoms 2 and 3" },
178 { "-rmin", FALSE0, etREAL, {&rmin}, "Minimum distance (nm)" },
179 { "-rmax", FALSE0, etREAL, {&rmax}, "Maximum distance (nm)" },
180 { "-cbin", FALSE0, etREAL, {&binwidth}, "Binwidth for the cosine" },
181 { "-rbin", FALSE0, etREAL, {&rbinw}, "Binwidth for r (nm)" },
182 { "-pbc", FALSE0, etBOOL, {&bPBC}, "Check PBC for the center of mass calculation. Only necessary when your reference group consists of several molecules." }
183 };
184
185 t_filenm fnm[] = {
186 { efTRX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
187 { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
188 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
189 { efXVG, NULL((void*)0), "sori.xvg", ffWRITE1<<2 },
190 { efXVG, "-no", "snor.xvg", ffWRITE1<<2 },
191 { efXVG, "-ro", "sord.xvg", ffWRITE1<<2 },
192 { efXVG, "-co", "scum.xvg", ffWRITE1<<2 },
193 { efXVG, "-rc", "scount.xvg", ffWRITE1<<2 }
194 };
195#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
196
197 if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_CAN_VIEW(1<<5) | PCA_BE_NICE(1<<13),
198 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))
199 {
200 return 0;
201 }
202
203 two_pi = 2/M_PI3.14159265358979323846;
204
205 bTPS = (opt2bSet("-s", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || !opt2bSet("-n", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm) || bCom);
206 if (bTPS)
207 {
208 read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &xtop, NULL((void*)0), box,
209 bCom);
210 }
211
212 /* get index groups */
213 printf("Select a group of reference particles and a solvent group:\n");
214 snew(grpname, 2)(grpname) = save_calloc("grpname", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 214, (2), sizeof(*(grpname)));
215 snew(index, 2)(index) = save_calloc("index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 215, (2), sizeof(*(index)));
216 snew(isize, 2)(isize) = save_calloc("isize", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 216, (2), sizeof(*(isize)));
217 if (bTPS)
218 {
219 get_index(&top.atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 2, isize, index, grpname);
220 }
221 else
222 {
223 get_index(NULL((void*)0), ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 2, isize, index, grpname);
224 }
225
226 if (bCom)
227 {
228 nrefgrp = 1;
229 nrefat = isize[0];
230 }
231 else
232 {
233 nrefgrp = isize[0];
234 nrefat = 1;
235 }
236
237 if (isize[1] % 3)
238 {
239 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 239, "The number of solvent atoms (%d) is not a multiple of 3",
240 isize[1]);
241 }
242
243 /* initialize reading trajectory: */
244 natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box);
245
246 rmin2 = sqr(rmin);
247 rmax2 = sqr(rmax);
248 rcut = 0.99*sqrt(max_cutoff2(guess_ePBC(box), box));
249 if (rcut == 0)
250 {
251 rcut = 10*rmax;
252 }
253 rcut2 = sqr(rcut);
254
255 invbw = 1/binwidth;
256 nbin1 = 1+(int)(2*invbw + 0.5);
257 nbin2 = 1+(int)(invbw + 0.5);
258
259 invrbw = 1/rbinw;
260
261 snew(hist1, nbin1)(hist1) = save_calloc("hist1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 261, (nbin1), sizeof(*(hist1)));
262 snew(hist2, nbin2)(hist2) = save_calloc("hist2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 262, (nbin2), sizeof(*(hist2)));
263 nrbin = 1+(int)(rcut/rbinw);
264 if (nrbin == 0)
265 {
266 nrbin = 1;
267 }
268 snew(histi1, nrbin)(histi1) = save_calloc("histi1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 268, (nrbin), sizeof(*(histi1)));
269 snew(histi2, nrbin)(histi2) = save_calloc("histi2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 269, (nrbin), sizeof(*(histi2)));
270 snew(histn, nrbin)(histn) = save_calloc("histn", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 270, (nrbin), sizeof(*(histn)));
271
272 ntot = 0;
273 nf = 0;
274 sum1 = 0;
275 sum2 = 0;
276
277 if (bTPS)
278 {
279 /* make molecules whole again */
280 gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
281 }
282 /* start analysis of trajectory */
283 do
284 {
285 if (bTPS)
286 {
287 /* make molecules whole again */
288 gmx_rmpbc(gpbc, natoms, box, x);
289 }
290
291 set_pbc(&pbc, ePBC, box);
292 n = 0;
293 inp = 0;
294 outp = 0;
Value stored to 'outp' is never read
295 for (p = 0; (p < nrefgrp); p++)
296 {
297 if (bCom)
298 {
299 calc_com_pbc(nrefat, &top, x, &pbc, index[0], xref, bPBC);
300 }
301 else
302 {
303 copy_rvec(x[index[0][p]], xref);
304 }
305
306 for (m = 0; m < isize[1]; m += 3)
307 {
308 sa0 = index[1][m];
309 sa1 = index[1][m+1];
310 sa2 = index[1][m+2];
311 range_check(sa0, 0, natoms)_range_check(sa0, 0, natoms, ((void*)0),"sa0", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 311);
312 range_check(sa1, 0, natoms)_range_check(sa1, 0, natoms, ((void*)0),"sa1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 312);
313 range_check(sa2, 0, natoms)_range_check(sa2, 0, natoms, ((void*)0),"sa2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 313);
314 pbc_dx(&pbc, x[sa0], xref, dx);
315 r2 = norm2(dx);
316 if (r2 < rcut2)
317 {
318 r = sqrt(r2);
319 if (!bVec23)
320 {
321 /* Determine the normal to the plain */
322 rvec_sub(x[sa1], x[sa0], dxh1);
323 rvec_sub(x[sa2], x[sa0], dxh2);
324 rvec_inc(dxh1, dxh2);
325 svmul(1/r, dx, dx);
326 unitv(dxh1, dxh1);
327 inp = iprod(dx, dxh1);
328 cprod(dxh1, dxh2, outer);
329 unitv(outer, outer);
330 outp = iprod(dx, outer);
331 }
332 else
333 {
334 /* Use the vector between the 2nd and 3rd atom */
335 rvec_sub(x[sa2], x[sa1], dxh2);
336 unitv(dxh2, dxh2);
337 outp = iprod(dx, dxh2)/r;
338 }
339 {
340 int ii = (int)(invrbw*r);
341 range_check(ii, 0, nrbin)_range_check(ii, 0, nrbin, ((void*)0),"ii", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 341);
342 histi1[ii] += inp;
343 histi2[ii] += 3*sqr(outp) - 1;
344 histn[ii]++;
345 }
346 if ((r2 >= rmin2) && (r2 < rmax2))
347 {
348 int ii1 = (int)(invbw*(inp + 1));
349 int ii2 = (int)(invbw*fabs(outp));
350
351 range_check(ii1, 0, nbin1)_range_check(ii1, 0, nbin1, ((void*)0),"ii1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 351);
352 range_check(ii2, 0, nbin2)_range_check(ii2, 0, nbin2, ((void*)0),"ii2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 352);
353 hist1[ii1]++;
354 hist2[ii2]++;
355 sum1 += inp;
356 sum2 += outp;
357 n++;
358 }
359 }
360 }
361 }
362 ntot += n;
363 nf++;
364
365 }
366 while (read_next_x(oenv, status, &t, x, box));
367
368 /* clean up */
369 sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_sorient.c"
, 369, (x));
370 close_trj(status);
371 gmx_rmpbc_done(gpbc);
372
373 /* Add the bin for the exact maximum to the previous bin */
374 hist1[nbin1-1] += hist1[nbin1];
375 hist2[nbin2-1] += hist2[nbin2];
376
377 nav = (real)ntot/(nrefgrp*nf);
378 normfac = invbw/ntot;
379
380 fprintf(stderrstderr, "Average nr of molecules between %g and %g nm: %.1f\n",
381 rmin, rmax, nav);
382 if (ntot > 0)
383 {
384 sum1 /= ntot;
385 sum2 /= ntot;
386 fprintf(stderrstderr, "Average cos(theta1) between %g and %g nm: %6.3f\n",
387 rmin, rmax, sum1);
388 fprintf(stderrstderr, "Average 3cos2(theta2)-1 between %g and %g nm: %6.3f\n",
389 rmin, rmax, sum2);
390 }
391
392 sprintf(str, "Solvent orientation between %g and %g nm", rmin, rmax);
393 fp = xvgropen(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "cos(\\8q\\4\\s1\\N)", "", oenv);
394 if (output_env_get_print_xvgr_codes(oenv))
395 {
396 fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
397 }
398 for (i = 0; i < nbin1; i++)
399 {
400 fprintf(fp, "%g %g\n", (i+0.5)*binwidth-1, 2*normfac*hist1[i]);
401 }
402 gmx_ffclose(fp);
403
404 sprintf(str, "Solvent normal orientation between %g and %g nm", rmin, rmax);
405 fp = xvgropen(opt2fn("-no", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "cos(\\8q\\4\\s2\\N)", "", oenv);
406 if (output_env_get_print_xvgr_codes(oenv))
407 {
408 fprintf(fp, "@ subtitle \"average shell size %.1f molecules\"\n", nav);
409 }
410 for (i = 0; i < nbin2; i++)
411 {
412 fprintf(fp, "%g %g\n", (i+0.5)*binwidth, normfac*hist2[i]);
413 }
414 gmx_ffclose(fp);
415
416
417 sprintf(str, "Solvent orientation");
418 fp = xvgropen(opt2fn("-ro", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "r (nm)", "", oenv);
419 if (output_env_get_print_xvgr_codes(oenv))
420 {
421 fprintf(fp, "@ subtitle \"as a function of distance\"\n");
422 }
423 xvgr_legend(fp, 2, legr, oenv);
424 for (i = 0; i < nrbin; i++)
425 {
426 fprintf(fp, "%g %g %g\n", (i+0.5)*rbinw,
427 histn[i] ? histi1[i]/histn[i] : 0,
428 histn[i] ? histi2[i]/histn[i] : 0);
429 }
430 gmx_ffclose(fp);
431
432 sprintf(str, "Cumulative solvent orientation");
433 fp = xvgropen(opt2fn("-co", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "r (nm)", "", oenv);
434 if (output_env_get_print_xvgr_codes(oenv))
435 {
436 fprintf(fp, "@ subtitle \"as a function of distance\"\n");
437 }
438 xvgr_legend(fp, 2, legc, oenv);
439 normfac = 1.0/(nrefgrp*nf);
440 c1 = 0;
441 c2 = 0;
442 fprintf(fp, "%g %g %g\n", 0.0, c1, c2);
443 for (i = 0; i < nrbin; i++)
444 {
445 c1 += histi1[i]*normfac;
446 c2 += histi2[i]*normfac;
447 fprintf(fp, "%g %g %g\n", (i+1)*rbinw, c1, c2);
448 }
449 gmx_ffclose(fp);
450
451 sprintf(str, "Solvent distribution");
452 fp = xvgropen(opt2fn("-rc", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), str, "r (nm)", "molecules/nm", oenv);
453 if (output_env_get_print_xvgr_codes(oenv))
454 {
455 fprintf(fp, "@ subtitle \"as a function of distance\"\n");
456 }
457 normfac = 1.0/(rbinw*nf);
458 for (i = 0; i < nrbin; i++)
459 {
460 fprintf(fp, "%g %g\n", (i+0.5)*rbinw, histn[i]*normfac);
461 }
462 gmx_ffclose(fp);
463
464 do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
465 do_view(oenv, opt2fn("-no", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0));
466 do_view(oenv, opt2fn("-ro", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
467 do_view(oenv, opt2fn("-co", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy");
468
469 return 0;
470}