File: | gromacs/gmxana/gmx_rotmat.c |
Location: | line 259, column 9 |
Description: | Array access results in a null pointer dereference |
1 | /* | |||
2 | * This file is part of the GROMACS molecular simulation package. | |||
3 | * | |||
4 | * Copyright (c) 2009,2010,2011,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 | #ifdef HAVE_CONFIG_H1 | |||
36 | #include <config.h> | |||
37 | #endif | |||
38 | ||||
39 | #include <math.h> | |||
40 | #include <string.h> | |||
41 | ||||
42 | #include "gromacs/commandline/pargs.h" | |||
43 | #include "typedefs.h" | |||
44 | #include "gromacs/utility/smalloc.h" | |||
45 | #include "macros.h" | |||
46 | #include "gromacs/math/vec.h" | |||
47 | #include "pbc.h" | |||
48 | #include "gromacs/utility/futil.h" | |||
49 | #include "index.h" | |||
50 | #include "mshift.h" | |||
51 | #include "gromacs/fileio/xvgr.h" | |||
52 | #include "viewit.h" | |||
53 | #include "rmpbc.h" | |||
54 | #include "gromacs/fileio/tpxio.h" | |||
55 | #include "gromacs/fileio/trxio.h" | |||
56 | #include "gmx_ana.h" | |||
57 | ||||
58 | #include "gromacs/math/do_fit.h" | |||
59 | ||||
60 | static void get_refx(output_env_t oenv, const char *trxfn, int nfitdim, int skip, | |||
61 | int gnx, int *index, | |||
62 | gmx_bool bMW, t_topology *top, int ePBC, rvec *x_ref) | |||
63 | { | |||
64 | int natoms, nfr_all, nfr, i, j, a, r, c, min_fr; | |||
65 | t_trxstatus *status; | |||
66 | real *ti, min_t; | |||
67 | double tot_mass, msd, *srmsd, min_srmsd, srmsd_tot; | |||
68 | rvec *x, **xi; | |||
69 | real xf; | |||
70 | matrix box, R; | |||
71 | real *w_rls; | |||
72 | gmx_rmpbc_t gpbc = NULL((void*)0); | |||
73 | ||||
74 | ||||
75 | nfr_all = 0; | |||
76 | nfr = 0; | |||
77 | snew(ti, 100)(ti) = save_calloc("ti", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 77, (100), sizeof(*(ti))); | |||
78 | snew(xi, 100)(xi) = save_calloc("xi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 78, (100), sizeof(*(xi))); | |||
79 | natoms = read_first_x(oenv, &status, trxfn, &ti[nfr], &x, box); | |||
80 | ||||
81 | snew(w_rls, gnx)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 81, (gnx), sizeof(*(w_rls))); | |||
82 | tot_mass = 0; | |||
83 | for (a = 0; a < gnx; a++) | |||
84 | { | |||
85 | if (index[a] >= natoms) | |||
86 | { | |||
87 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 87, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[a]+1, natoms); | |||
88 | } | |||
89 | w_rls[a] = (bMW ? top->atoms.atom[index[a]].m : 1.0); | |||
90 | tot_mass += w_rls[a]; | |||
91 | } | |||
92 | gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms); | |||
93 | ||||
94 | do | |||
95 | { | |||
96 | if (nfr_all % skip == 0) | |||
97 | { | |||
98 | gmx_rmpbc(gpbc, natoms, box, x); | |||
99 | snew(xi[nfr], gnx)(xi[nfr]) = save_calloc("xi[nfr]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 99, (gnx), sizeof(*(xi[nfr]))); | |||
100 | for (i = 0; i < gnx; i++) | |||
101 | { | |||
102 | copy_rvec(x[index[i]], xi[nfr][i]); | |||
103 | } | |||
104 | reset_x(gnx, NULL((void*)0), gnx, NULL((void*)0), xi[nfr], w_rls); | |||
105 | nfr++; | |||
106 | if (nfr % 100 == 0) | |||
107 | { | |||
108 | srenew(ti, nfr+100)(ti) = save_realloc("ti", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 108, (ti), (nfr+100), sizeof(*(ti))); | |||
109 | srenew(xi, nfr+100)(xi) = save_realloc("xi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 109, (xi), (nfr+100), sizeof(*(xi))); | |||
110 | } | |||
111 | } | |||
112 | nfr_all++; | |||
113 | } | |||
114 | while (read_next_x(oenv, status, &ti[nfr], x, box)); | |||
115 | close_trj(status); | |||
116 | sfree(x)save_free("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 116, (x)); | |||
117 | ||||
118 | gmx_rmpbc_done(gpbc); | |||
119 | ||||
120 | snew(srmsd, nfr)(srmsd) = save_calloc("srmsd", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 120, (nfr), sizeof(*(srmsd))); | |||
121 | for (i = 0; i < nfr; i++) | |||
122 | { | |||
123 | printf("\rProcessing frame %d of %d", i, nfr); | |||
124 | for (j = i+1; j < nfr; j++) | |||
125 | { | |||
126 | calc_fit_R(nfitdim, gnx, w_rls, xi[i], xi[j], R); | |||
127 | ||||
128 | msd = 0; | |||
129 | for (a = 0; a < gnx; a++) | |||
130 | { | |||
131 | for (r = 0; r < DIM3; r++) | |||
132 | { | |||
133 | xf = 0; | |||
134 | for (c = 0; c < DIM3; c++) | |||
135 | { | |||
136 | xf += R[r][c]*xi[j][a][c]; | |||
137 | } | |||
138 | msd += w_rls[a]*sqr(xi[i][a][r] - xf); | |||
139 | } | |||
140 | } | |||
141 | msd /= tot_mass; | |||
142 | srmsd[i] += sqrt(msd); | |||
143 | srmsd[j] += sqrt(msd); | |||
144 | } | |||
145 | sfree(xi[i])save_free("xi[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 145, (xi[i])); | |||
146 | } | |||
147 | printf("\n"); | |||
148 | sfree(w_rls)save_free("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 148, (w_rls)); | |||
149 | ||||
150 | min_srmsd = GMX_REAL_MAX3.40282346E+38; | |||
151 | min_fr = -1; | |||
152 | min_t = -1; | |||
153 | srmsd_tot = 0; | |||
154 | for (i = 0; i < nfr; i++) | |||
155 | { | |||
156 | srmsd[i] /= (nfr - 1); | |||
157 | if (srmsd[i] < min_srmsd) | |||
158 | { | |||
159 | min_srmsd = srmsd[i]; | |||
160 | min_fr = i; | |||
161 | min_t = ti[i]; | |||
162 | } | |||
163 | srmsd_tot += srmsd[i]; | |||
164 | } | |||
165 | sfree(srmsd)save_free("srmsd", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 165, (srmsd)); | |||
166 | ||||
167 | printf("Average RMSD between all structures: %.3f\n", srmsd_tot/nfr); | |||
168 | printf("Structure with lowest RMSD to all others: time %g, av. RMSD %.3f\n", | |||
169 | min_t, min_srmsd); | |||
170 | ||||
171 | for (a = 0; a < gnx; a++) | |||
172 | { | |||
173 | copy_rvec(xi[min_fr][a], x_ref[index[a]]); | |||
174 | } | |||
175 | ||||
176 | sfree(xi)save_free("xi", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 176, (xi)); | |||
177 | } | |||
178 | ||||
179 | int gmx_rotmat(int argc, char *argv[]) | |||
180 | { | |||
181 | const char *desc[] = { | |||
182 | "[THISMODULE] plots the rotation matrix required for least squares fitting", | |||
183 | "a conformation onto the reference conformation provided with", | |||
184 | "[TT]-s[tt]. Translation is removed before fitting.", | |||
185 | "The output are the three vectors that give the new directions", | |||
186 | "of the x, y and z directions of the reference conformation,", | |||
187 | "for example: (zx,zy,zz) is the orientation of the reference", | |||
188 | "z-axis in the trajectory frame.", | |||
189 | "[PAR]", | |||
190 | "This tool is useful for, for instance,", | |||
191 | "determining the orientation of a molecule", | |||
192 | "at an interface, possibly on a trajectory produced with", | |||
193 | "[TT]gmx trjconv -fit rotxy+transxy[tt] to remove the rotation", | |||
194 | "in the [IT]x-y[it] plane.", | |||
195 | "[PAR]", | |||
196 | "Option [TT]-ref[tt] determines a reference structure for fitting,", | |||
197 | "instead of using the structure from [TT]-s[tt]. The structure with", | |||
198 | "the lowest sum of RMSD's to all other structures is used.", | |||
199 | "Since the computational cost of this procedure grows with", | |||
200 | "the square of the number of frames, the [TT]-skip[tt] option", | |||
201 | "can be useful. A full fit or only a fit in the [IT]x-y[it] plane can", | |||
202 | "be performed.", | |||
203 | "[PAR]", | |||
204 | "Option [TT]-fitxy[tt] fits in the [IT]x-y[it] plane before determining", | |||
205 | "the rotation matrix." | |||
206 | }; | |||
207 | const char *reffit[] = | |||
208 | { NULL((void*)0), "none", "xyz", "xy", NULL((void*)0) }; | |||
209 | static int skip = 1; | |||
210 | static gmx_bool bFitXY = FALSE0, bMW = TRUE1; | |||
211 | t_pargs pa[] = { | |||
212 | { "-ref", FALSE0, etENUM, {reffit}, | |||
213 | "Determine the optimal reference structure" }, | |||
214 | { "-skip", FALSE0, etINT, {&skip}, | |||
215 | "Use every nr-th frame for [TT]-ref[tt]" }, | |||
216 | { "-fitxy", FALSE0, etBOOL, {&bFitXY}, | |||
217 | "Fit the x/y rotation before determining the rotation" }, | |||
218 | { "-mw", FALSE0, etBOOL, {&bMW}, | |||
219 | "Use mass weighted fitting" } | |||
220 | }; | |||
221 | FILE *out; | |||
222 | t_trxstatus *status; | |||
223 | t_topology top; | |||
224 | int ePBC; | |||
225 | rvec *x_ref, *x; | |||
226 | matrix box, R; | |||
227 | real t; | |||
228 | int natoms, i; | |||
229 | char *grpname, title[256]; | |||
230 | int gnx; | |||
231 | gmx_rmpbc_t gpbc = NULL((void*)0); | |||
232 | atom_id *index; | |||
233 | output_env_t oenv; | |||
234 | real *w_rls; | |||
235 | const char *leg[] = { "xx", "xy", "xz", "yx", "yy", "yz", "zx", "zy", "zz" }; | |||
236 | #define NLEG((int)(sizeof(leg)/sizeof((leg)[0]))) asize(leg)((int)(sizeof(leg)/sizeof((leg)[0]))) | |||
237 | t_filenm fnm[] = { | |||
238 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, | |||
239 | { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, | |||
240 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, | |||
241 | { efXVG, NULL((void*)0), "rotmat", ffWRITE1<<2 } | |||
242 | }; | |||
243 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) | |||
244 | ||||
245 | 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), | |||
| ||||
246 | 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)) | |||
247 | { | |||
248 | return 0; | |||
249 | } | |||
250 | ||||
251 | read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), title, &top, &ePBC, &x_ref, NULL((void*)0), box, bMW); | |||
252 | ||||
253 | gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr); | |||
254 | ||||
255 | gmx_rmpbc(gpbc, top.atoms.nr, box, x_ref); | |||
256 | ||||
257 | get_index(&top.atoms, ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &gnx, &index, &grpname); | |||
258 | ||||
259 | if (reffit[0][0] != 'n') | |||
| ||||
260 | { | |||
261 | get_refx(oenv, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), reffit[0][2] == 'z' ? 3 : 2, skip, | |||
262 | gnx, index, bMW, &top, ePBC, x_ref); | |||
263 | } | |||
264 | ||||
265 | natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box); | |||
266 | ||||
267 | snew(w_rls, natoms)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 267, (natoms), sizeof(*(w_rls))); | |||
268 | for (i = 0; i < gnx; i++) | |||
269 | { | |||
270 | if (index[i] >= natoms) | |||
271 | { | |||
272 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rotmat.c" , 272, "Atom index (%d) is larger than the number of atoms in the trajecory (%d)", index[i]+1, natoms); | |||
273 | } | |||
274 | w_rls[index[i]] = (bMW ? top.atoms.atom[index[i]].m : 1.0); | |||
275 | } | |||
276 | ||||
277 | if (reffit[0][0] == 'n') | |||
278 | { | |||
279 | reset_x(gnx, index, natoms, NULL((void*)0), x_ref, w_rls); | |||
280 | } | |||
281 | ||||
282 | out = xvgropen(ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), | |||
283 | "Fit matrix", "Time (ps)", "", oenv); | |||
284 | xvgr_legend(out, NLEG((int)(sizeof(leg)/sizeof((leg)[0]))), leg, oenv); | |||
285 | ||||
286 | do | |||
287 | { | |||
288 | gmx_rmpbc(gpbc, natoms, box, x); | |||
289 | ||||
290 | reset_x(gnx, index, natoms, NULL((void*)0), x, w_rls); | |||
291 | ||||
292 | if (bFitXY) | |||
293 | { | |||
294 | do_fit_ndim(2, natoms, w_rls, x_ref, x); | |||
295 | } | |||
296 | ||||
297 | calc_fit_R(DIM3, natoms, w_rls, x_ref, x, R); | |||
298 | ||||
299 | fprintf(out, | |||
300 | "%7g %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f\n", | |||
301 | t, | |||
302 | R[XX0][XX0], R[XX0][YY1], R[XX0][ZZ2], | |||
303 | R[YY1][XX0], R[YY1][YY1], R[YY1][ZZ2], | |||
304 | R[ZZ2][XX0], R[ZZ2][YY1], R[ZZ2][ZZ2]); | |||
305 | } | |||
306 | while (read_next_x(oenv, status, &t, x, box)); | |||
307 | ||||
308 | gmx_rmpbc_done(gpbc); | |||
309 | ||||
310 | close_trj(status); | |||
311 | ||||
312 | gmx_ffclose(out); | |||
313 | ||||
314 | do_view(oenv, ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-nxy"); | |||
315 | ||||
316 | return 0; | |||
317 | } |