File: | gromacs/gmxana/gmx_rotacf.c |
Location: | line 161, column 5 |
Description: | Value stored to 't1' is never read |
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 <string.h> |
43 | |
44 | #include "physics.h" |
45 | #include "typedefs.h" |
46 | #include "gromacs/utility/smalloc.h" |
47 | #include "gromacs/utility/futil.h" |
48 | #include "gromacs/commandline/pargs.h" |
49 | #include "index.h" |
50 | #include "macros.h" |
51 | #include "gromacs/utility/fatalerror.h" |
52 | #include "gstat.h" |
53 | #include "gromacs/math/vec.h" |
54 | #include "viewit.h" |
55 | #include "gmx_ana.h" |
56 | #include "gromacs/fileio/trxio.h" |
57 | |
58 | |
59 | int gmx_rotacf(int argc, char *argv[]) |
60 | { |
61 | const char *desc[] = { |
62 | "[THISMODULE] calculates the rotational correlation function", |
63 | "for molecules. Atom triplets (i,j,k) must be given in the index", |
64 | "file, defining two vectors ij and jk. The rotational ACF", |
65 | "is calculated as the autocorrelation function of the vector", |
66 | "n = ij x jk, i.e. the cross product of the two vectors.", |
67 | "Since three atoms span a plane, the order of the three atoms", |
68 | "does not matter. Optionally, by invoking the [TT]-d[tt] switch, you can", |
69 | "calculate the rotational correlation function for linear molecules", |
70 | "by specifying atom pairs (i,j) in the index file.", |
71 | "[PAR]", |
72 | "EXAMPLES[PAR]", |
73 | "[TT]gmx rotacf -P 1 -nparm 2 -fft -n index -o rotacf-x-P1", |
74 | "-fa expfit-x-P1 -beginfit 2.5 -endfit 20.0[tt][PAR]", |
75 | "This will calculate the rotational correlation function using a first", |
76 | "order Legendre polynomial of the angle of a vector defined by the index", |
77 | "file. The correlation function will be fitted from 2.5 ps until 20.0 ps", |
78 | "to a two-parameter exponential." |
79 | }; |
80 | static gmx_bool bVec = FALSE0, bAver = TRUE1; |
81 | |
82 | t_pargs pa[] = { |
83 | { "-d", FALSE0, etBOOL, {&bVec}, |
84 | "Use index doublets (vectors) for correlation function instead of triplets (planes)" }, |
85 | { "-aver", FALSE0, etBOOL, {&bAver}, |
86 | "Average over molecules" } |
87 | }; |
88 | |
89 | t_trxstatus *status; |
90 | int isize; |
91 | atom_id *index; |
92 | char *grpname; |
93 | rvec *x, *x_s; |
94 | matrix box; |
95 | real **c1; |
96 | rvec xij, xjk, n; |
97 | int i, m, teller, n_alloc, natoms, nvec, ai, aj, ak; |
98 | unsigned long mode; |
99 | real t, t0, t1, dt; |
100 | gmx_rmpbc_t gpbc = NULL((void*)0); |
101 | t_topology *top; |
102 | int ePBC; |
103 | t_filenm fnm[] = { |
104 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, |
105 | { efTPX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
106 | { efNDX, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
107 | { efXVG, "-o", "rotacf", ffWRITE1<<2 } |
108 | }; |
109 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
110 | int npargs; |
111 | t_pargs *ppa; |
112 | |
113 | output_env_t oenv; |
114 | |
115 | npargs = asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))); |
116 | ppa = add_acf_pargs(&npargs, pa); |
117 | |
118 | 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), |
119 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, npargs, ppa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, 0, NULL((void*)0), &oenv)) |
120 | { |
121 | return 0; |
122 | } |
123 | |
124 | rd_index(ftp2fn(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &isize, &index, &grpname); |
125 | |
126 | if (bVec) |
127 | { |
128 | nvec = isize/2; |
129 | } |
130 | else |
131 | { |
132 | nvec = isize/3; |
133 | } |
134 | |
135 | if (((isize % 3) != 0) && !bVec) |
136 | { |
137 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotacf.c" , 137, "number of index elements not multiple of 3, " |
138 | "these can not be atom triplets\n"); |
139 | } |
140 | if (((isize % 2) != 0) && bVec) |
141 | { |
142 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotacf.c" , 142, "number of index elements not multiple of 2, " |
143 | "these can not be atom doublets\n"); |
144 | } |
145 | |
146 | top = read_top(ftp2fn(efTPX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &ePBC); |
147 | |
148 | snew(c1, nvec)(c1) = save_calloc("c1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotacf.c" , 148, (nvec), sizeof(*(c1))); |
149 | for (i = 0; (i < nvec); i++) |
150 | { |
151 | c1[i] = NULL((void*)0); |
152 | } |
153 | n_alloc = 0; |
154 | |
155 | natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box); |
156 | snew(x_s, natoms)(x_s) = save_calloc("x_s", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotacf.c" , 156, (natoms), sizeof(*(x_s))); |
157 | |
158 | gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms); |
159 | |
160 | /* Start the loop over frames */ |
161 | t1 = t0 = t; |
Value stored to 't1' is never read | |
162 | teller = 0; |
163 | do |
164 | { |
165 | if (teller >= n_alloc) |
166 | { |
167 | n_alloc += 100; |
168 | for (i = 0; (i < nvec); i++) |
169 | { |
170 | srenew(c1[i], DIM*n_alloc)(c1[i]) = save_realloc("c1[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotacf.c" , 170, (c1[i]), (3*n_alloc), sizeof(*(c1[i]))); |
171 | } |
172 | } |
173 | t1 = t; |
174 | |
175 | /* Remove periodicity */ |
176 | gmx_rmpbc_copy(gpbc, natoms, box, x, x_s); |
177 | |
178 | /* Compute crossproducts for all vectors, if triplets. |
179 | * else, just get the vectors in case of doublets. |
180 | */ |
181 | if (bVec == FALSE0) |
182 | { |
183 | for (i = 0; (i < nvec); i++) |
184 | { |
185 | ai = index[3*i]; |
186 | aj = index[3*i+1]; |
187 | ak = index[3*i+2]; |
188 | rvec_sub(x_s[ai], x_s[aj], xij); |
189 | rvec_sub(x_s[aj], x_s[ak], xjk); |
190 | cprod(xij, xjk, n); |
191 | for (m = 0; (m < DIM3); m++) |
192 | { |
193 | c1[i][DIM3*teller+m] = n[m]; |
194 | } |
195 | } |
196 | } |
197 | else |
198 | { |
199 | for (i = 0; (i < nvec); i++) |
200 | { |
201 | ai = index[2*i]; |
202 | aj = index[2*i+1]; |
203 | rvec_sub(x_s[ai], x_s[aj], n); |
204 | for (m = 0; (m < DIM3); m++) |
205 | { |
206 | c1[i][DIM3*teller+m] = n[m]; |
207 | } |
208 | } |
209 | } |
210 | /* Increment loop counter */ |
211 | teller++; |
212 | } |
213 | while (read_next_x(oenv, status, &t, x, box)); |
214 | close_trj(status); |
215 | fprintf(stderrstderr, "\nDone with trajectory\n"); |
216 | |
217 | gmx_rmpbc_done(gpbc); |
218 | |
219 | |
220 | /* Autocorrelation function */ |
221 | if (teller < 2) |
222 | { |
223 | fprintf(stderrstderr, "Not enough frames for correlation function\n"); |
224 | } |
225 | else |
226 | { |
227 | dt = (t1 - t0)/(teller-1); |
228 | |
229 | mode = eacVector(1<<2); |
230 | |
231 | do_autocorr(ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), oenv, "Rotational Correlation Function", |
232 | teller, nvec, c1, dt, mode, bAver); |
233 | } |
234 | |
235 | do_view(oenv, ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
236 | |
237 | return 0; |
238 | } |