Bug Summary

File:gromacs/gmxana/gmx_dyecoupl.c
Location:line 137, column 5
Description:Value stored to 'in_ndxfile' is never read

Annotated Source Code

1/*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2012,2013,2014, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
8 *
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
13 *
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
18 *
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
23 *
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
31 *
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
34 */
35#include "copyrite.h"
36#include "gromacs/fileio/filenm.h"
37#include "macros.h"
38#include "pbc.h"
39#include "gromacs/utility/smalloc.h"
40#include "gromacs/commandline/pargs.h"
41#include "gromacs/math/vec.h"
42#include "gromacs/fileio/xvgr.h"
43#include "gromacs/fileio/trxio.h"
44
45#include "gromacs/utility/fatalerror.h"
46
47int gmx_dyecoupl(int argc, char *argv[])
48{
49 const char *desc[] =
50 {
51 "[THISMODULE] extracts dye dynamics from trajectory files.",
52 "Currently, R and kappa^2 between dyes is extracted for (F)RET",
53 "simulations with assumed dipolar coupling as in the Foerster equation.",
54 "It further allows the calculation of R(t) and kappa^2(t), R and",
55 "kappa^2 histograms and averages, as well as the instantaneous FRET",
56 "efficiency E(t) for a specified Foerster radius R_0 (switch [TT]-R0[tt]).",
57 "The input dyes have to be whole (see res and mol pbc options",
58 "in [TT]trjconv[tt]).",
59 "The dye transition dipole moment has to be defined by at least",
60 "a single atom pair, however multiple atom pairs can be provided ",
61 "in the index file. The distance R is calculated on the basis of",
62 "the COMs of the given atom pairs.",
63 "The [TT]-pbcdist[tt] option calculates distances to the nearest periodic",
64 "image instead to the distance in the box. This works however only,"
65 "for periodic boundaries in all 3 dimensions.",
66 "The [TT]-norm[tt] option (area-) normalizes the histograms."
67 };
68
69 static gmx_bool bPBCdist = FALSE0, bNormHist = FALSE0;
70 int histbins = 50;
71 output_env_t oenv;
72 real R0 = -1;
73
74 t_pargs pa[] =
75 {
76 { "-pbcdist", FALSE0, etBOOL, { &bPBCdist }, "Distance R based on PBC" },
77 { "-norm", FALSE0, etBOOL, { &bNormHist }, "Normalize histograms" },
78 { "-bins", FALSE0, etINT, {&histbins}, "# of histogram bins" },
79 { "-R0", FALSE0, etREAL, {&R0}, "Foerster radius including kappa^2=2/3 in nm" }
80 };
81#define NPA((int)(sizeof(pa)/sizeof((pa)[0]))) asize(pa)((int)(sizeof(pa)/sizeof((pa)[0])))
82
83 t_filenm fnm[] =
84 {
85 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
86 { efNDX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
87 { efXVG, "-ot", "rkappa", ffOPTWR(1<<2| 1<<3) },
88 { efXVG, "-oe", "insteff", ffOPTWR(1<<2| 1<<3) },
89 { efDAT, "-o", "rkappa", ffOPTWR(1<<2| 1<<3) },
90 { efXVG, "-rhist", "rhist", ffOPTWR(1<<2| 1<<3) },
91 { efXVG, "-khist", "khist", ffOPTWR(1<<2| 1<<3) }
92 };
93#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
94
95
96 const char *in_trajfile, *in_ndxfile, *out_xvgrkfile = NULL((void*)0), *out_xvginstefffile = NULL((void*)0), *out_xvgrhistfile = NULL((void*)0), *out_xvgkhistfile = NULL((void*)0), *out_datfile = NULL((void*)0);
97 gmx_bool bHaveFirstFrame, bHaveNextFrame, indexOK = TRUE1;
98 int ndon, nacc;
99 atom_id *donindex, *accindex;
100 char *grpnm;
101 t_atoms *atoms = NULL((void*)0);
102 t_trxstatus *status;
103 t_trxframe fr;
104
105 int flags;
106 int allocblock = 1000;
107 real histexpand = 1e-6;
108 rvec donvec, accvec, donpos, accpos, dist, distnorm;
109 int natoms;
110
111 /*we rely on PBC autodetection (...currently)*/
112 int ePBC = -1;
113
114 real *rvalues = NULL((void*)0), *kappa2values = NULL((void*)0), *rhist = NULL((void*)0), *khist = NULL((void*)0);
115 t_pbc *pbc = NULL((void*)0);
116 int i, bin;
117 FILE *rkfp = NULL((void*)0), *rhfp = NULL((void*)0), *khfp = NULL((void*)0), *datfp = NULL((void*)0), *iefp = NULL((void*)0);
118 gmx_bool bRKout, bRhistout, bKhistout, bDatout, bInstEffout, grident;
119
120 const char *rkleg[2] = { "R", "\\f{Symbol}k\\f{}\\S2\\N" };
121 const char *rhleg[1] = { "p(R)" };
122 const char *khleg[1] = { "p(\\f{Symbol}k\\f{}\\S2\\N)" };
123 const char *ieleg[1] = { "E\\sRET\\N(t)" };
124
125 real R, kappa2, insteff, Rs = 0., kappa2s = 0., insteffs = 0., rmax, rmin, kmin = 0., kmax = 4.,
126 rrange, krange, rincr, kincr, Rfrac;
127 int rkcount = 0, rblocksallocated = 0, kblocksallocated = 0;
128
129 if (!parse_common_args(&argc, argv, PCA_CAN_BEGIN(1<<6) | PCA_CAN_END(1<<7) | PCA_CAN_VIEW(1<<5) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13), 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, NULL((void*)0), &oenv))
130 {
131 return 0;
132 }
133
134
135 /* Check command line options for filenames and set bool flags when switch used*/
136 in_trajfile = opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
137 in_ndxfile = opt2fn("-n", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
Value stored to 'in_ndxfile' is never read
138 out_xvgrkfile = opt2fn("-ot", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
139 out_xvgrhistfile = opt2fn("-rhist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
140 out_xvgkhistfile = opt2fn("-khist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
141 out_xvginstefffile = opt2fn("-oe", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
142 out_datfile = opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
143
144 bRKout = opt2bSet("-ot", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
145 bRhistout = opt2bSet("-rhist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
146 bKhistout = opt2bSet("-khist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
147 bDatout = opt2bSet("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
148 bInstEffout = opt2bSet("-oe", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
149
150
151 /* PBC warning. */
152 if (bPBCdist)
153 {
154 printf("Calculating distances to periodic image.\n");
155 printf("Be careful! This produces only valid results for PBC in all three dimensions\n");
156 }
157
158
159 if (bInstEffout && R0 <= 0.)
160 {
161 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 161, "You have to specify R0 and R0 has to be larger than 0 nm.\n\n");
162 }
163
164 printf("Select group with donor atom pairs defining the transition moment\n");
165 get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &ndon, &donindex, &grpnm);
166
167 printf("Select group with acceptor atom pairs defining the transition moment\n");
168 get_index(atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &nacc, &accindex, &grpnm);
169
170 /*check if groups are identical*/
171 grident = TRUE1;
172
173 if (ndon == nacc)
174 {
175 for (i = 0; i < nacc; i++)
176 {
177 if (accindex[i] != donindex[i])
178 {
179 grident = FALSE0;
180 break;
181 }
182 }
183 }
184
185 if (grident)
186 {
187 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 187, "Donor and acceptor group are identical. This makes no sense.");
188 }
189
190 printf("Reading first frame\n");
191 /* open trx file for reading */
192 flags = 0;
193 flags = flags | TRX_READ_X(1<<0);
194 bHaveFirstFrame = read_first_frame(oenv, &status, in_trajfile, &fr, flags);
195
196 if (bHaveFirstFrame)
197 {
198 printf("First frame is OK\n");
199 natoms = fr.natoms;
200 if ((ndon % 2 != 0) || (nacc % 2 != 0))
201 {
202 indexOK = FALSE0;
203 }
204 else
205 {
206 for (i = 0; i < ndon; i++)
207 {
208 if (donindex[i] >= natoms)
209 {
210 indexOK = FALSE0;
211 }
212 }
213 for (i = 0; i < nacc; i++)
214 {
215 if (accindex[i] >= natoms)
216 {
217 indexOK = FALSE0;
218 }
219 }
220 }
221
222 if (indexOK)
223 {
224
225 if (bDatout)
226 {
227 datfp = fopen(out_datfile, "w");
228 }
229
230 if (bRKout)
231 {
232 rkfp = xvgropen(out_xvgrkfile,
233 "Distance and \\f{Symbol}k\\f{}\\S2\\N trajectory",
234 "Time (ps)", "Distance (nm) / \\f{Symbol}k\\f{}\\S2\\N",
235 oenv);
236 xvgr_legend(rkfp, 2, rkleg, oenv);
237 }
238
239 if (bInstEffout)
240 {
241 iefp = xvgropen(out_xvginstefffile,
242 "Instantaneous RET Efficiency",
243 "Time (ps)", "RET Efficiency",
244 oenv);
245 xvgr_legend(iefp, 1, ieleg, oenv);
246 }
247
248
249 if (bRhistout)
250 {
251 snew(rvalues, allocblock)(rvalues) = save_calloc("rvalues", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 251, (allocblock), sizeof(*(rvalues)));
252 rblocksallocated += 1;
253 snew(rhist, histbins)(rhist) = save_calloc("rhist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 253, (histbins), sizeof(*(rhist)));
254 }
255
256 if (bKhistout)
257 {
258 snew(kappa2values, allocblock)(kappa2values) = save_calloc("kappa2values", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 258, (allocblock), sizeof(*(kappa2values)));
259 kblocksallocated += 1;
260 snew(khist, histbins)(khist) = save_calloc("khist", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 260, (histbins), sizeof(*(khist)));
261 }
262
263 do
264 {
265 clear_rvec(donvec);
266 clear_rvec(accvec);
267 clear_rvec(donpos);
268 clear_rvec(accpos);
269 for (i = 0; i < ndon / 2; i++)
270 {
271 rvec_sub(donvec, fr.x[donindex[2 * i]], donvec);
272 rvec_add(donvec, fr.x[donindex[2 * i + 1]], donvec);
273 rvec_add(donpos, fr.x[donindex[2 * i]], donpos);
274 rvec_add(donpos, fr.x[donindex[2 * i + 1]], donpos);
275 }
276
277 for (i = 0; i < nacc / 2; i++)
278 {
279 rvec_sub(accvec, fr.x[accindex[2 * i]], accvec);
280 rvec_add(accvec, fr.x[accindex[2 * i + 1]], accvec);
281 rvec_add(accpos, fr.x[accindex[2 * i]], accpos);
282 rvec_add(accpos, fr.x[accindex[2 * i + 1]], accpos);
283 }
284
285 unitv(donvec, donvec);
286 unitv(accvec, accvec);
287
288 svmul((real) 1. / ndon, donpos, donpos);
289 svmul((real) 1. / nacc, accpos, accpos);
290
291 if (bPBCdist)
292 {
293 set_pbc(pbc, ePBC, fr.box);
294 pbc_dx(pbc, donpos, accpos, dist);
295 }
296 else
297 {
298 rvec_sub(donpos, accpos, dist);
299 }
300
301 unitv(dist, distnorm);
302 R = norm(dist);
303 kappa2 = iprod(donvec, accvec)- 3.* (iprod(donvec, distnorm) * iprod(distnorm, accvec));
304 kappa2 *= kappa2;
305 if (R0 > 0)
306 {
307 Rfrac = R/R0;
308 insteff = 1/(1+(Rfrac*Rfrac*Rfrac*Rfrac*Rfrac*Rfrac)*2/3/kappa2);
309 insteffs += insteff;
310
311 if (bInstEffout)
312 {
313 fprintf(iefp, "%12.7f %12.7f\n", fr.time, insteff);
314 }
315 }
316
317
318 Rs += R;
319 kappa2s += kappa2;
320 rkcount++;
321
322 if (bRKout)
323 {
324 fprintf(rkfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
325 }
326
327 if (bDatout)
328 {
329 fprintf(datfp, "%12.7f %12.7f %12.7f\n", fr.time, R, kappa2);
330 }
331
332 if (bRhistout)
333 {
334 rvalues[rkcount-1] = R;
335 if (rkcount % allocblock == 0)
336 {
337 srenew(rvalues, allocblock*(rblocksallocated+1))(rvalues) = save_realloc("rvalues", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 337, (rvalues), (allocblock*(rblocksallocated+1)), sizeof(*
(rvalues)));
338 rblocksallocated += 1;
339 }
340 }
341
342 if (bKhistout)
343 {
344 kappa2values[rkcount-1] = kappa2;
345 if (rkcount % allocblock == 0)
346 {
347 srenew(kappa2values, allocblock*(kblocksallocated+1))(kappa2values) = save_realloc("kappa2values", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 347, (kappa2values), (allocblock*(kblocksallocated+1)), sizeof
(*(kappa2values)));
348 kblocksallocated += 1;
349 }
350 }
351
352 bHaveNextFrame = read_next_frame(oenv, status, &fr);
353 }
354 while (bHaveNextFrame);
355
356 if (bRKout)
357 {
358 gmx_ffclose(rkfp);
359 }
360
361 if (bDatout)
362 {
363 gmx_ffclose(datfp);
364 }
365
366 if (bInstEffout)
367 {
368 gmx_ffclose(iefp);
369 }
370
371
372 if (bRhistout)
373 {
374 printf("Writing R-Histogram\n");
375 rmin = rvalues[0];
376 rmax = rvalues[0];
377 for (i = 1; i < rkcount; i++)
378 {
379 if (rvalues[i] < rmin)
380 {
381 rmin = rvalues[i];
382 }
383 else if (rvalues[i] > rmax)
384 {
385 rmax = rvalues[i];
386 }
387 }
388 rmin -= histexpand;
389 rmax += histexpand;
390
391 rrange = rmax - rmin;
392 rincr = rrange / histbins;
393
394 for (i = 1; i < rkcount; i++)
395 {
396 bin = (int) ((rvalues[i] - rmin) / rincr);
397 rhist[bin] += 1;
398 }
399 if (bNormHist)
400 {
401 for (i = 0; i < histbins; i++)
402 {
403 rhist[i] /= rkcount * rrange/histbins;
404 }
405 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
406 "R (nm)", "Normalized Probability", oenv);
407 }
408 else
409 {
410 rhfp = xvgropen(out_xvgrhistfile, "Distance Distribution",
411 "R (nm)", "Probability", oenv);
412 }
413 xvgr_legend(rhfp, 1, rhleg, oenv);
414 for (i = 0; i < histbins; i++)
415 {
416 fprintf(rhfp, "%12.7f %12.7f\n", (i + 0.5) * rincr + rmin,
417 rhist[i]);
418 }
419 gmx_ffclose(rhfp);
420 }
421
422 if (bKhistout)
423 {
424 printf("Writing kappa^2-Histogram\n");
425 krange = kmax - kmin;
426 kincr = krange / histbins;
427
428 for (i = 1; i < rkcount; i++)
429 {
430 bin = (int) ((kappa2values[i] - kmin) / kincr);
431 khist[bin] += 1;
432 }
433 if (bNormHist)
434 {
435 for (i = 0; i < histbins; i++)
436 {
437 khist[i] /= rkcount * krange/histbins;
438 }
439 khfp = xvgropen(out_xvgkhistfile,
440 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
441 "\\f{Symbol}k\\f{}\\S2\\N",
442 "Normalized Probability", oenv);
443 }
444 else
445 {
446 khfp = xvgropen(out_xvgkhistfile,
447 "\\f{Symbol}k\\f{}\\S2\\N Distribution",
448 "\\f{Symbol}k\\f{}\\S2\\N", "Probability", oenv);
449 }
450 xvgr_legend(khfp, 1, khleg, oenv);
451 for (i = 0; i < histbins; i++)
452 {
453 fprintf(khfp, "%12.7f %12.7f\n", (i + 0.5) * kincr + kmin,
454 khist[i]);
455 }
456 gmx_ffclose(khfp);
457 }
458
459 printf("\nAverages:\n");
460 printf("R_avg = %8.4f nm\nKappa^2 = %8.4f\n", Rs / rkcount,
461 kappa2s / rkcount);
462 if (R0 > 0)
463 {
464 printf("E_RETavg = %8.4f\n", insteffs / rkcount);
465 }
466 please_cite(stdoutstdout, "Hoefling2011");
467 }
468 else
469 {
470 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 470, "Index file invalid, check your index file for correct pairs.\n");
471 }
472 }
473 else
474 {
475 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_dyecoupl.c"
, 475, "Could not read first frame of the trajectory.\n");
476 }
477
478 return 0;
479}