File: | gromacs/gmxana/gmx_rms.c |
Location: | line 564, column 17 |
Description: | Value stored to 'iatom' 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 <stdlib.h> |
43 | |
44 | #include "gromacs/utility/smalloc.h" |
45 | #include "macros.h" |
46 | #include "typedefs.h" |
47 | #include "gromacs/fileio/xvgr.h" |
48 | #include "copyrite.h" |
49 | #include "gromacs/commandline/pargs.h" |
50 | #include "gromacs/math/vec.h" |
51 | #include "index.h" |
52 | #include "gromacs/utility/fatalerror.h" |
53 | #include "gromacs/utility/futil.h" |
54 | #include "princ.h" |
55 | #include "rmpbc.h" |
56 | #include "gromacs/fileio/matio.h" |
57 | #include "gromacs/fileio/tpxio.h" |
58 | #include "gromacs/fileio/trxio.h" |
59 | #include "cmat.h" |
60 | #include "viewit.h" |
61 | #include "gmx_ana.h" |
62 | |
63 | #include "gromacs/math/do_fit.h" |
64 | |
65 | static void norm_princ(t_atoms *atoms, int isize, atom_id *index, int natoms, |
66 | rvec *x) |
67 | { |
68 | int i, m; |
69 | rvec princ, vec; |
70 | |
71 | /* equalize principal components: */ |
72 | /* orient principal axes, get principal components */ |
73 | orient_princ(atoms, isize, index, natoms, x, NULL((void*)0), princ); |
74 | |
75 | /* calc our own principal components */ |
76 | clear_rvec(vec); |
77 | for (m = 0; m < DIM3; m++) |
78 | { |
79 | for (i = 0; i < isize; i++) |
80 | { |
81 | vec[m] += sqr(x[index[i]][m]); |
82 | } |
83 | vec[m] = sqrt(vec[m] / isize); |
84 | /* calculate scaling constants */ |
85 | vec[m] = 1 / (sqrt(3) * vec[m]); |
86 | } |
87 | |
88 | /* scale coordinates */ |
89 | for (i = 0; i < natoms; i++) |
90 | { |
91 | for (m = 0; m < DIM3; m++) |
92 | { |
93 | x[i][m] *= vec[m]; |
94 | } |
95 | } |
96 | } |
97 | |
98 | int gmx_rms(int argc, char *argv[]) |
99 | { |
100 | const char *desc[] = |
101 | { |
102 | "[THISMODULE] compares two structures by computing the root mean square", |
103 | "deviation (RMSD), the size-independent [GRK]rho[grk] similarity parameter", |
104 | "([TT]rho[tt]) or the scaled [GRK]rho[grk] ([TT]rhosc[tt]), ", |
105 | "see Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).", |
106 | "This is selected by [TT]-what[tt].[PAR]" |
107 | |
108 | "Each structure from a trajectory ([TT]-f[tt]) is compared to a", |
109 | "reference structure. The reference structure", |
110 | "is taken from the structure file ([TT]-s[tt]).[PAR]", |
111 | |
112 | "With option [TT]-mir[tt] also a comparison with the mirror image of", |
113 | "the reference structure is calculated.", |
114 | "This is useful as a reference for 'significant' values, see", |
115 | "Maiorov & Crippen, Proteins [BB]22[bb], 273 (1995).[PAR]", |
116 | |
117 | "Option [TT]-prev[tt] produces the comparison with a previous frame", |
118 | "the specified number of frames ago.[PAR]", |
119 | |
120 | "Option [TT]-m[tt] produces a matrix in [TT].xpm[tt] format of", |
121 | "comparison values of each structure in the trajectory with respect to", |
122 | "each other structure. This file can be visualized with for instance", |
123 | "[TT]xv[tt] and can be converted to postscript with [gmx-xpm2ps].[PAR]", |
124 | |
125 | "Option [TT]-fit[tt] controls the least-squares fitting of", |
126 | "the structures on top of each other: complete fit (rotation and", |
127 | "translation), translation only, or no fitting at all.[PAR]", |
128 | |
129 | "Option [TT]-mw[tt] controls whether mass weighting is done or not.", |
130 | "If you select the option (default) and ", |
131 | "supply a valid [TT].tpr[tt] file masses will be taken from there, ", |
132 | "otherwise the masses will be deduced from the [TT]atommass.dat[tt] file in", |
133 | "[TT]GMXLIB[tt]. This is fine for proteins, but not", |
134 | "necessarily for other molecules. A default mass of 12.011 amu (carbon)", |
135 | "is assigned to unknown atoms. You can check whether this happend by", |
136 | "turning on the [TT]-debug[tt] flag and inspecting the log file.[PAR]", |
137 | |
138 | "With [TT]-f2[tt], the 'other structures' are taken from a second", |
139 | "trajectory, this generates a comparison matrix of one trajectory", |
140 | "versus the other.[PAR]", |
141 | |
142 | "Option [TT]-bin[tt] does a binary dump of the comparison matrix.[PAR]", |
143 | |
144 | "Option [TT]-bm[tt] produces a matrix of average bond angle deviations", |
145 | "analogously to the [TT]-m[tt] option. Only bonds between atoms in the", |
146 | "comparison group are considered." |
147 | }; |
148 | static gmx_bool bPBC = TRUE1, bFitAll = TRUE1, bSplit = FALSE0; |
149 | static gmx_bool bDeltaLog = FALSE0; |
150 | static int prev = 0, freq = 1, freq2 = 1, nlevels = 80, avl = 0; |
151 | static real rmsd_user_max = -1, rmsd_user_min = -1, bond_user_max = -1, |
152 | bond_user_min = -1, delta_maxy = 0.0; |
153 | /* strings and things for selecting difference method */ |
154 | enum |
155 | { |
156 | ewSel, ewRMSD, ewRho, ewRhoSc, ewNR |
157 | }; |
158 | int ewhat; |
159 | const char *what[ewNR + 1] = |
160 | { NULL((void*)0), "rmsd", "rho", "rhosc", NULL((void*)0) }; |
161 | const char *whatname[ewNR] = |
162 | { NULL((void*)0), "RMSD", "Rho", "Rho sc" }; |
163 | const char *whatlabel[ewNR] = |
164 | { NULL((void*)0), "RMSD (nm)", "Rho", "Rho sc" }; |
165 | const char *whatxvgname[ewNR] = |
166 | { NULL((void*)0), "RMSD", "\\8r\\4", "\\8r\\4\\ssc\\N" }; |
167 | const char *whatxvglabel[ewNR] = |
168 | { NULL((void*)0), "RMSD (nm)", "\\8r\\4", "\\8r\\4\\ssc\\N" }; |
169 | /* strings and things for fitting methods */ |
170 | enum |
171 | { |
172 | efSel, efFit, efReset, efNone, efNR |
173 | }; |
174 | int efit; |
175 | const char *fit[efNR + 1] = |
176 | { NULL((void*)0), "rot+trans", "translation", "none", NULL((void*)0) }; |
177 | const char *fitgraphlabel[efNR + 1] = |
178 | { NULL((void*)0), "lsq fit", "translational fit", "no fit" }; |
179 | static int nrms = 1; |
180 | static gmx_bool bMassWeighted = TRUE1; |
181 | t_pargs pa[] = |
182 | { |
183 | { "-what", FALSE0, etENUM, |
184 | { what }, "Structural difference measure" }, |
185 | { "-pbc", FALSE0, etBOOL, |
186 | { &bPBC }, "PBC check" }, |
187 | { "-fit", FALSE0, etENUM, |
188 | { fit }, "Fit to reference structure" }, |
189 | { "-prev", FALSE0, etINT, |
190 | { &prev }, "Compare with previous frame" }, |
191 | { "-split", FALSE0, etBOOL, |
192 | { &bSplit }, "Split graph where time is zero" }, |
193 | { "-fitall", FALSE0, etBOOL, |
194 | { &bFitAll }, "HIDDENFit all pairs of structures in matrix" }, |
195 | { "-skip", FALSE0, etINT, |
196 | { &freq }, "Only write every nr-th frame to matrix" }, |
197 | { "-skip2", FALSE0, etINT, |
198 | { &freq2 }, "Only write every nr-th frame to matrix" }, |
199 | { "-max", FALSE0, etREAL, |
200 | { &rmsd_user_max }, "Maximum level in comparison matrix" }, |
201 | { "-min", FALSE0, etREAL, |
202 | { &rmsd_user_min }, "Minimum level in comparison matrix" }, |
203 | { "-bmax", FALSE0, etREAL, |
204 | { &bond_user_max }, "Maximum level in bond angle matrix" }, |
205 | { "-bmin", FALSE0, etREAL, |
206 | { &bond_user_min }, "Minimum level in bond angle matrix" }, |
207 | { "-mw", FALSE0, etBOOL, |
208 | { &bMassWeighted }, "Use mass weighting for superposition" }, |
209 | { "-nlevels", FALSE0, etINT, |
210 | { &nlevels }, "Number of levels in the matrices" }, |
211 | { "-ng", FALSE0, etINT, |
212 | { &nrms }, "Number of groups to compute RMS between" }, |
213 | { "-dlog", FALSE0, etBOOL, |
214 | { &bDeltaLog }, |
215 | "HIDDENUse a log x-axis in the delta t matrix" }, |
216 | { "-dmax", FALSE0, etREAL, |
217 | { &delta_maxy }, "HIDDENMaximum level in delta matrix" }, |
218 | { "-aver", FALSE0, etINT, |
219 | { &avl }, |
220 | "HIDDENAverage over this distance in the RMSD matrix" } |
221 | }; |
222 | int natoms_trx, natoms_trx2, natoms; |
223 | int i, j, k, m, teller, teller2, tel_mat, tel_mat2; |
224 | #define NFRAME5000 5000 |
225 | int maxframe = NFRAME5000, maxframe2 = NFRAME5000; |
226 | real t, *w_rls, *w_rms, *w_rls_m = NULL((void*)0), *w_rms_m = NULL((void*)0); |
227 | gmx_bool bNorm, bAv, bFreq2, bFile2, bMat, bBond, bDelta, bMirror, bMass; |
228 | gmx_bool bFit, bReset; |
229 | t_topology top; |
230 | int ePBC; |
231 | t_iatom *iatom = NULL((void*)0); |
232 | |
233 | matrix box; |
234 | rvec *x, *xp, *xm = NULL((void*)0), **mat_x = NULL((void*)0), **mat_x2, *mat_x2_j = NULL((void*)0), vec1, |
235 | vec2; |
236 | t_trxstatus *status; |
237 | char buf[256], buf2[256]; |
238 | int ncons = 0; |
239 | FILE *fp; |
240 | real rlstot = 0, **rls, **rlsm = NULL((void*)0), *time, *time2, *rlsnorm = NULL((void*)0), |
241 | **rmsd_mat = NULL((void*)0), **bond_mat = NULL((void*)0), *axis, *axis2, *del_xaxis, |
242 | *del_yaxis, rmsd_max, rmsd_min, rmsd_avg, bond_max, bond_min, ang; |
243 | real **rmsdav_mat = NULL((void*)0), av_tot, weight, weight_tot; |
244 | real **delta = NULL((void*)0), delta_max, delta_scalex = 0, delta_scaley = 0, |
245 | *delta_tot; |
246 | int delta_xsize = 0, del_lev = 100, mx, my, abs_my; |
247 | gmx_bool bA1, bA2, bPrev, bTop, *bInMat = NULL((void*)0); |
248 | int ifit, *irms, ibond = 0, *ind_bond1 = NULL((void*)0), *ind_bond2 = NULL((void*)0), n_ind_m = |
249 | 0; |
250 | atom_id *ind_fit, **ind_rms, *ind_m = NULL((void*)0), *rev_ind_m = NULL((void*)0), *ind_rms_m = |
251 | NULL((void*)0); |
252 | char *gn_fit, **gn_rms; |
253 | t_rgb rlo, rhi; |
254 | output_env_t oenv; |
255 | gmx_rmpbc_t gpbc = NULL((void*)0); |
256 | |
257 | t_filenm fnm[] = |
258 | { |
259 | { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 }, |
260 | { efTRX, "-f", NULL((void*)0), ffREAD1<<1 }, |
261 | { efTRX, "-f2", NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
262 | { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) }, |
263 | { efXVG, NULL((void*)0), "rmsd", ffWRITE1<<2 }, |
264 | { efXVG, "-mir", "rmsdmir", ffOPTWR(1<<2| 1<<3) }, |
265 | { efXVG, "-a", "avgrp", ffOPTWR(1<<2| 1<<3) }, |
266 | { efXVG, "-dist", "rmsd-dist", ffOPTWR(1<<2| 1<<3) }, |
267 | { efXPM, "-m", "rmsd", ffOPTWR(1<<2| 1<<3) }, |
268 | { efDAT, "-bin", "rmsd", ffOPTWR(1<<2| 1<<3) }, |
269 | { efXPM, "-bm", "bond", ffOPTWR(1<<2| 1<<3) } |
270 | }; |
271 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
272 | |
273 | if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_CAN_VIEW(1<<5) |
274 | | PCA_BE_NICE(1<<13), 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), |
275 | &oenv)) |
276 | { |
277 | return 0; |
278 | } |
279 | /* parse enumerated options: */ |
280 | ewhat = nenum(what); |
281 | if (ewhat == ewRho || ewhat == ewRhoSc) |
282 | { |
283 | please_cite(stdoutstdout, "Maiorov95"); |
284 | } |
285 | efit = nenum(fit); |
286 | bFit = efit == efFit; |
287 | bReset = efit == efReset; |
288 | if (bFit) |
289 | { |
290 | bReset = TRUE1; /* for fit, reset *must* be set */ |
291 | } |
292 | else |
293 | { |
294 | bFitAll = FALSE0; |
295 | } |
296 | |
297 | /* mark active cmdline options */ |
298 | bMirror = opt2bSet("-mir", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); /* calc RMSD vs mirror of ref. */ |
299 | bFile2 = opt2bSet("-f2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
300 | bMat = opt2bSet("-m", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
301 | bBond = opt2bSet("-bm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
302 | bDelta = (delta_maxy > 0); /* calculate rmsd vs delta t matrix from * |
303 | * your RMSD matrix (hidden option */ |
304 | bNorm = opt2bSet("-a", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
305 | bFreq2 = opt2parg_bSet("-skip2", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa); |
306 | if (freq <= 0) |
307 | { |
308 | fprintf(stderrstderr, "The number of frames to skip is <= 0. " |
309 | "Writing out all frames.\n\n"); |
310 | freq = 1; |
311 | } |
312 | if (!bFreq2) |
313 | { |
314 | freq2 = freq; |
315 | } |
316 | else if (bFile2 && freq2 <= 0) |
317 | { |
318 | fprintf(stderrstderr, |
319 | "The number of frames to skip in second trajectory is <= 0.\n" |
320 | " Writing out all frames.\n\n"); |
321 | freq2 = 1; |
322 | } |
323 | |
324 | bPrev = (prev > 0); |
325 | if (bPrev) |
326 | { |
327 | prev = abs(prev); |
328 | if (freq != 1) |
329 | { |
330 | fprintf(stderrstderr, "WARNING: option -skip also applies to -prev\n"); |
331 | } |
332 | } |
333 | |
334 | if (bFile2 && !bMat && !bBond) |
335 | { |
336 | fprintf( |
337 | stderrstderr, |
338 | "WARNING: second trajectory (-f2) useless when not calculating matrix (-m/-bm),\n" |
339 | " will not read from %s\n", opt2fn("-f2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), |
340 | fnm)); |
341 | bFile2 = FALSE0; |
342 | } |
343 | |
344 | if (bDelta) |
345 | { |
346 | bMat = TRUE1; |
347 | if (bFile2) |
348 | { |
349 | fprintf(stderrstderr, |
350 | "WARNING: second trajectory (-f2) useless when making delta matrix,\n" |
351 | " will not read from %s\n", opt2fn("-f2", |
352 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)); |
353 | bFile2 = FALSE0; |
354 | } |
355 | } |
356 | |
357 | bTop = read_tps_conf(ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), buf, &top, &ePBC, &xp, |
358 | NULL((void*)0), box, TRUE1); |
359 | snew(w_rls, top.atoms.nr)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 359, (top.atoms.nr), sizeof(*(w_rls))); |
360 | snew(w_rms, top.atoms.nr)(w_rms) = save_calloc("w_rms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 360, (top.atoms.nr), sizeof(*(w_rms))); |
361 | |
362 | if (!bTop && bBond) |
363 | { |
364 | fprintf(stderrstderr, |
365 | "WARNING: Need a run input file for bond angle matrix,\n" |
366 | " will not calculate bond angle matrix.\n"); |
367 | bBond = FALSE0; |
368 | } |
369 | |
370 | if (bReset) |
371 | { |
372 | fprintf(stderrstderr, "Select group for %s fit\n", bFit ? "least squares" |
373 | : "translational"); |
374 | get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), 1, &ifit, |
375 | &ind_fit, &gn_fit); |
376 | } |
377 | else |
378 | { |
379 | ifit = 0; |
380 | } |
381 | |
382 | if (bReset) |
383 | { |
384 | if (bFit && ifit < 3) |
385 | { |
386 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 386, "Need >= 3 points to fit!\n" ); |
387 | } |
388 | |
389 | bMass = FALSE0; |
390 | for (i = 0; i < ifit; i++) |
391 | { |
392 | if (bMassWeighted) |
393 | { |
394 | w_rls[ind_fit[i]] = top.atoms.atom[ind_fit[i]].m; |
395 | } |
396 | else |
397 | { |
398 | w_rls[ind_fit[i]] = 1; |
399 | } |
400 | bMass = bMass || (top.atoms.atom[ind_fit[i]].m != 0); |
401 | } |
402 | if (!bMass) |
403 | { |
404 | fprintf(stderrstderr, "All masses in the fit group are 0, using masses of 1\n"); |
405 | for (i = 0; i < ifit; i++) |
406 | { |
407 | w_rls[ind_fit[i]] = 1; |
408 | } |
409 | } |
410 | } |
411 | |
412 | if (bMat || bBond) |
413 | { |
414 | nrms = 1; |
415 | } |
416 | |
417 | snew(gn_rms, nrms)(gn_rms) = save_calloc("gn_rms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 417, (nrms), sizeof(*(gn_rms))); |
418 | snew(ind_rms, nrms)(ind_rms) = save_calloc("ind_rms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 418, (nrms), sizeof(*(ind_rms))); |
419 | snew(irms, nrms)(irms) = save_calloc("irms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 419, (nrms), sizeof(*(irms))); |
420 | |
421 | fprintf(stderrstderr, "Select group%s for %s calculation\n", |
422 | (nrms > 1) ? "s" : "", whatname[ewhat]); |
423 | get_index(&(top.atoms), ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), |
424 | nrms, irms, ind_rms, gn_rms); |
425 | |
426 | if (bNorm) |
427 | { |
428 | snew(rlsnorm, irms[0])(rlsnorm) = save_calloc("rlsnorm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 428, (irms[0]), sizeof(*(rlsnorm))); |
429 | } |
430 | snew(rls, nrms)(rls) = save_calloc("rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 430, (nrms), sizeof(*(rls))); |
431 | for (j = 0; j < nrms; j++) |
432 | { |
433 | snew(rls[j], maxframe)(rls[j]) = save_calloc("rls[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 433, (maxframe), sizeof(*(rls[j]))); |
434 | } |
435 | if (bMirror) |
436 | { |
437 | snew(rlsm, nrms)(rlsm) = save_calloc("rlsm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 437, (nrms), sizeof(*(rlsm))); |
438 | for (j = 0; j < nrms; j++) |
439 | { |
440 | snew(rlsm[j], maxframe)(rlsm[j]) = save_calloc("rlsm[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 440, (maxframe), sizeof(*(rlsm[j]))); |
441 | } |
442 | } |
443 | snew(time, maxframe)(time) = save_calloc("time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 443, (maxframe), sizeof(*(time))); |
444 | for (j = 0; j < nrms; j++) |
445 | { |
446 | bMass = FALSE0; |
447 | for (i = 0; i < irms[j]; i++) |
448 | { |
449 | if (bMassWeighted) |
450 | { |
451 | w_rms[ind_rms[j][i]] = top.atoms.atom[ind_rms[j][i]].m; |
452 | } |
453 | else |
454 | { |
455 | w_rms[ind_rms[j][i]] = 1.0; |
456 | } |
457 | bMass = bMass || (top.atoms.atom[ind_rms[j][i]].m != 0); |
458 | } |
459 | if (!bMass) |
460 | { |
461 | fprintf(stderrstderr, "All masses in group %d are 0, using masses of 1\n", j); |
462 | for (i = 0; i < irms[j]; i++) |
463 | { |
464 | w_rms[ind_rms[j][i]] = 1; |
465 | } |
466 | } |
467 | } |
468 | /* Prepare reference frame */ |
469 | if (bPBC) |
470 | { |
471 | gpbc = gmx_rmpbc_init(&top.idef, ePBC, top.atoms.nr); |
472 | gmx_rmpbc(gpbc, top.atoms.nr, box, xp); |
473 | } |
474 | if (bReset) |
475 | { |
476 | reset_x(ifit, ind_fit, top.atoms.nr, NULL((void*)0), xp, w_rls); |
477 | } |
478 | if (bMirror) |
479 | { |
480 | /* generate reference structure mirror image: */ |
481 | snew(xm, top.atoms.nr)(xm) = save_calloc("xm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 481, (top.atoms.nr), sizeof(*(xm))); |
482 | for (i = 0; i < top.atoms.nr; i++) |
483 | { |
484 | copy_rvec(xp[i], xm[i]); |
485 | xm[i][XX0] = -xm[i][XX0]; |
486 | } |
487 | } |
488 | if (ewhat == ewRhoSc) |
489 | { |
490 | norm_princ(&top.atoms, ifit, ind_fit, top.atoms.nr, xp); |
491 | } |
492 | |
493 | /* read first frame */ |
494 | natoms_trx = read_first_x(oenv, &status, opt2fn("-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box); |
495 | if (natoms_trx != top.atoms.nr) |
496 | { |
497 | fprintf(stderrstderr, |
498 | "\nWARNING: topology has %d atoms, whereas trajectory has %d\n", |
499 | top.atoms.nr, natoms_trx); |
500 | } |
501 | natoms = min(top.atoms.nr, natoms_trx)(((top.atoms.nr) < (natoms_trx)) ? (top.atoms.nr) : (natoms_trx ) ); |
502 | if (bMat || bBond || bPrev) |
503 | { |
504 | snew(mat_x, NFRAME)(mat_x) = save_calloc("mat_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 504, (5000), sizeof(*(mat_x))); |
505 | |
506 | if (bPrev) |
507 | { |
508 | /* With -prev we use all atoms for simplicity */ |
509 | n_ind_m = natoms; |
510 | } |
511 | else |
512 | { |
513 | /* Check which atoms we need (fit/rms) */ |
514 | snew(bInMat, natoms)(bInMat) = save_calloc("bInMat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 514, (natoms), sizeof(*(bInMat))); |
515 | for (i = 0; i < ifit; i++) |
516 | { |
517 | bInMat[ind_fit[i]] = TRUE1; |
518 | } |
519 | n_ind_m = ifit; |
520 | for (i = 0; i < irms[0]; i++) |
521 | { |
522 | if (!bInMat[ind_rms[0][i]]) |
523 | { |
524 | bInMat[ind_rms[0][i]] = TRUE1; |
525 | n_ind_m++; |
526 | } |
527 | } |
528 | } |
529 | /* Make an index of needed atoms */ |
530 | snew(ind_m, n_ind_m)(ind_m) = save_calloc("ind_m", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 530, (n_ind_m), sizeof(*(ind_m))); |
531 | snew(rev_ind_m, natoms)(rev_ind_m) = save_calloc("rev_ind_m", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 531, (natoms), sizeof(*(rev_ind_m))); |
532 | j = 0; |
533 | for (i = 0; i < natoms; i++) |
534 | { |
535 | if (bPrev || bInMat[i]) |
536 | { |
537 | ind_m[j] = i; |
538 | rev_ind_m[i] = j; |
539 | j++; |
540 | } |
541 | } |
542 | snew(w_rls_m, n_ind_m)(w_rls_m) = save_calloc("w_rls_m", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 542, (n_ind_m), sizeof(*(w_rls_m))); |
543 | snew(ind_rms_m, irms[0])(ind_rms_m) = save_calloc("ind_rms_m", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 543, (irms[0]), sizeof(*(ind_rms_m))); |
544 | snew(w_rms_m, n_ind_m)(w_rms_m) = save_calloc("w_rms_m", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 544, (n_ind_m), sizeof(*(w_rms_m))); |
545 | for (i = 0; i < ifit; i++) |
546 | { |
547 | w_rls_m[rev_ind_m[ind_fit[i]]] = w_rls[ind_fit[i]]; |
548 | } |
549 | for (i = 0; i < irms[0]; i++) |
550 | { |
551 | ind_rms_m[i] = rev_ind_m[ind_rms[0][i]]; |
552 | w_rms_m[ind_rms_m[i]] = w_rms[ind_rms[0][i]]; |
553 | } |
554 | sfree(bInMat)save_free("bInMat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 554, (bInMat)); |
555 | } |
556 | |
557 | if (bBond) |
558 | { |
559 | ncons = 0; |
560 | for (k = 0; k < F_NRE; k++) |
561 | { |
562 | if (IS_CHEMBOND(k)(interaction_function[(k)].nratoms == 2 && (interaction_function [(k)].flags & 1<<3))) |
563 | { |
564 | iatom = top.idef.il[k].iatoms; |
Value stored to 'iatom' is never read | |
565 | ncons += top.idef.il[k].nr/3; |
566 | } |
567 | } |
568 | fprintf(stderrstderr, "Found %d bonds in topology\n", ncons); |
569 | snew(ind_bond1, ncons)(ind_bond1) = save_calloc("ind_bond1", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 569, (ncons), sizeof(*(ind_bond1))); |
570 | snew(ind_bond2, ncons)(ind_bond2) = save_calloc("ind_bond2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 570, (ncons), sizeof(*(ind_bond2))); |
571 | ibond = 0; |
572 | for (k = 0; k < F_NRE; k++) |
573 | { |
574 | if (IS_CHEMBOND(k)(interaction_function[(k)].nratoms == 2 && (interaction_function [(k)].flags & 1<<3))) |
575 | { |
576 | iatom = top.idef.il[k].iatoms; |
577 | ncons = top.idef.il[k].nr/3; |
578 | for (i = 0; i < ncons; i++) |
579 | { |
580 | bA1 = FALSE0; |
581 | bA2 = FALSE0; |
582 | for (j = 0; j < irms[0]; j++) |
583 | { |
584 | if (iatom[3*i+1] == ind_rms[0][j]) |
585 | { |
586 | bA1 = TRUE1; |
587 | } |
588 | if (iatom[3*i+2] == ind_rms[0][j]) |
589 | { |
590 | bA2 = TRUE1; |
591 | } |
592 | } |
593 | if (bA1 && bA2) |
594 | { |
595 | ind_bond1[ibond] = rev_ind_m[iatom[3*i+1]]; |
596 | ind_bond2[ibond] = rev_ind_m[iatom[3*i+2]]; |
597 | ibond++; |
598 | } |
599 | } |
600 | } |
601 | } |
602 | fprintf(stderrstderr, "Using %d bonds for bond angle matrix\n", ibond); |
603 | if (ibond == 0) |
604 | { |
605 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 605, "0 bonds found"); |
606 | } |
607 | } |
608 | |
609 | /* start looping over frames: */ |
610 | tel_mat = 0; |
611 | teller = 0; |
612 | do |
613 | { |
614 | if (bPBC) |
615 | { |
616 | gmx_rmpbc(gpbc, natoms, box, x); |
617 | } |
618 | |
619 | if (bReset) |
620 | { |
621 | reset_x(ifit, ind_fit, natoms, NULL((void*)0), x, w_rls); |
622 | } |
623 | if (ewhat == ewRhoSc) |
624 | { |
625 | norm_princ(&top.atoms, ifit, ind_fit, natoms, x); |
626 | } |
627 | |
628 | if (bFit) |
629 | { |
630 | /*do the least squares fit to original structure*/ |
631 | do_fit(natoms, w_rls, xp, x); |
632 | } |
633 | |
634 | if (teller % freq == 0) |
635 | { |
636 | /* keep frame for matrix calculation */ |
637 | if (bMat || bBond || bPrev) |
638 | { |
639 | if (tel_mat >= NFRAME5000) |
640 | { |
641 | srenew(mat_x, tel_mat+1)(mat_x) = save_realloc("mat_x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 641, (mat_x), (tel_mat+1), sizeof(*(mat_x))); |
642 | } |
643 | snew(mat_x[tel_mat], n_ind_m)(mat_x[tel_mat]) = save_calloc("mat_x[tel_mat]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 643, (n_ind_m), sizeof(*(mat_x[tel_mat]))); |
644 | for (i = 0; i < n_ind_m; i++) |
645 | { |
646 | copy_rvec(x[ind_m[i]], mat_x[tel_mat][i]); |
647 | } |
648 | } |
649 | tel_mat++; |
650 | } |
651 | |
652 | /*calculate energy of root_least_squares*/ |
653 | if (bPrev) |
654 | { |
655 | j = tel_mat-prev-1; |
656 | if (j < 0) |
657 | { |
658 | j = 0; |
659 | } |
660 | for (i = 0; i < n_ind_m; i++) |
661 | { |
662 | copy_rvec(mat_x[j][i], xp[ind_m[i]]); |
663 | } |
664 | if (bReset) |
665 | { |
666 | reset_x(ifit, ind_fit, natoms, NULL((void*)0), xp, w_rls); |
667 | } |
668 | if (bFit) |
669 | { |
670 | do_fit(natoms, w_rls, x, xp); |
671 | } |
672 | } |
673 | for (j = 0; (j < nrms); j++) |
674 | { |
675 | rls[j][teller] = |
676 | calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp); |
677 | } |
678 | if (bNorm) |
679 | { |
680 | for (j = 0; (j < irms[0]); j++) |
681 | { |
682 | rlsnorm[j] += |
683 | calc_similar_ind(ewhat != ewRMSD, 1, &(ind_rms[0][j]), w_rms, x, xp); |
684 | } |
685 | } |
686 | |
687 | if (bMirror) |
688 | { |
689 | if (bFit) |
690 | { |
691 | /*do the least squares fit to mirror of original structure*/ |
692 | do_fit(natoms, w_rls, xm, x); |
693 | } |
694 | |
695 | for (j = 0; j < nrms; j++) |
696 | { |
697 | rlsm[j][teller] = |
698 | calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xm); |
699 | } |
700 | } |
701 | time[teller] = output_env_conv_time(oenv, t); |
702 | |
703 | teller++; |
704 | if (teller >= maxframe) |
705 | { |
706 | maxframe += NFRAME5000; |
707 | srenew(time, maxframe)(time) = save_realloc("time", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 707, (time), (maxframe), sizeof(*(time))); |
708 | for (j = 0; (j < nrms); j++) |
709 | { |
710 | srenew(rls[j], maxframe)(rls[j]) = save_realloc("rls[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 710, (rls[j]), (maxframe), sizeof(*(rls[j]))); |
711 | } |
712 | if (bMirror) |
713 | { |
714 | for (j = 0; (j < nrms); j++) |
715 | { |
716 | srenew(rlsm[j], maxframe)(rlsm[j]) = save_realloc("rlsm[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 716, (rlsm[j]), (maxframe), sizeof(*(rlsm[j]))); |
717 | } |
718 | } |
719 | } |
720 | } |
721 | while (read_next_x(oenv, status, &t, x, box)); |
722 | close_trj(status); |
723 | |
724 | if (bFile2) |
725 | { |
726 | snew(time2, maxframe2)(time2) = save_calloc("time2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 726, (maxframe2), sizeof(*(time2))); |
727 | |
728 | fprintf(stderrstderr, "\nWill read second trajectory file\n"); |
729 | snew(mat_x2, NFRAME)(mat_x2) = save_calloc("mat_x2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 729, (5000), sizeof(*(mat_x2))); |
730 | natoms_trx2 = |
731 | read_first_x(oenv, &status, opt2fn("-f2", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), &t, &x, box); |
732 | if (natoms_trx2 != natoms_trx) |
733 | { |
734 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 734, |
735 | "Second trajectory (%d atoms) does not match the first one" |
736 | " (%d atoms)", natoms_trx2, natoms_trx); |
737 | } |
738 | tel_mat2 = 0; |
739 | teller2 = 0; |
740 | do |
741 | { |
742 | if (bPBC) |
743 | { |
744 | gmx_rmpbc(gpbc, natoms, box, x); |
745 | } |
746 | |
747 | if (bReset) |
748 | { |
749 | reset_x(ifit, ind_fit, natoms, NULL((void*)0), x, w_rls); |
750 | } |
751 | if (ewhat == ewRhoSc) |
752 | { |
753 | norm_princ(&top.atoms, ifit, ind_fit, natoms, x); |
754 | } |
755 | |
756 | if (bFit) |
757 | { |
758 | /*do the least squares fit to original structure*/ |
759 | do_fit(natoms, w_rls, xp, x); |
760 | } |
761 | |
762 | if (teller2 % freq2 == 0) |
763 | { |
764 | /* keep frame for matrix calculation */ |
765 | if (bMat) |
766 | { |
767 | if (tel_mat2 >= NFRAME5000) |
768 | { |
769 | srenew(mat_x2, tel_mat2+1)(mat_x2) = save_realloc("mat_x2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 769, (mat_x2), (tel_mat2+1), sizeof(*(mat_x2))); |
770 | } |
771 | snew(mat_x2[tel_mat2], n_ind_m)(mat_x2[tel_mat2]) = save_calloc("mat_x2[tel_mat2]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 771, (n_ind_m), sizeof(*(mat_x2[tel_mat2]))); |
772 | for (i = 0; i < n_ind_m; i++) |
773 | { |
774 | copy_rvec(x[ind_m[i]], mat_x2[tel_mat2][i]); |
775 | } |
776 | } |
777 | tel_mat2++; |
778 | } |
779 | |
780 | time2[teller2] = output_env_conv_time(oenv, t); |
781 | |
782 | teller2++; |
783 | if (teller2 >= maxframe2) |
784 | { |
785 | maxframe2 += NFRAME5000; |
786 | srenew(time2, maxframe2)(time2) = save_realloc("time2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 786, (time2), (maxframe2), sizeof(*(time2))); |
787 | } |
788 | } |
789 | while (read_next_x(oenv, status, &t, x, box)); |
790 | close_trj(status); |
791 | } |
792 | else |
793 | { |
794 | mat_x2 = mat_x; |
795 | time2 = time; |
796 | tel_mat2 = tel_mat; |
797 | freq2 = freq; |
798 | } |
799 | gmx_rmpbc_done(gpbc); |
800 | |
801 | if (bMat || bBond) |
802 | { |
803 | /* calculate RMS matrix */ |
804 | fprintf(stderrstderr, "\n"); |
805 | if (bMat) |
806 | { |
807 | fprintf(stderrstderr, "Building %s matrix, %dx%d elements\n", |
808 | whatname[ewhat], tel_mat, tel_mat2); |
809 | snew(rmsd_mat, tel_mat)(rmsd_mat) = save_calloc("rmsd_mat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 809, (tel_mat), sizeof(*(rmsd_mat))); |
810 | } |
811 | if (bBond) |
812 | { |
813 | fprintf(stderrstderr, "Building bond angle matrix, %dx%d elements\n", |
814 | tel_mat, tel_mat2); |
815 | snew(bond_mat, tel_mat)(bond_mat) = save_calloc("bond_mat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 815, (tel_mat), sizeof(*(bond_mat))); |
816 | } |
817 | snew(axis, tel_mat)(axis) = save_calloc("axis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 817, (tel_mat), sizeof(*(axis))); |
818 | snew(axis2, tel_mat2)(axis2) = save_calloc("axis2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 818, (tel_mat2), sizeof(*(axis2))); |
819 | rmsd_max = 0; |
820 | if (bFile2) |
821 | { |
822 | rmsd_min = 1e10; |
823 | } |
824 | else |
825 | { |
826 | rmsd_min = 0; |
827 | } |
828 | rmsd_avg = 0; |
829 | bond_max = 0; |
830 | bond_min = 1e10; |
831 | for (j = 0; j < tel_mat2; j++) |
832 | { |
833 | axis2[j] = time2[freq2*j]; |
834 | } |
835 | if (bDelta) |
836 | { |
837 | if (bDeltaLog) |
838 | { |
839 | delta_scalex = 8.0/log(2.0); |
840 | delta_xsize = (int)(log(tel_mat/2)*delta_scalex+0.5)+1; |
841 | } |
842 | else |
843 | { |
844 | delta_xsize = tel_mat/2; |
845 | } |
846 | delta_scaley = 1.0/delta_maxy; |
847 | snew(delta, delta_xsize)(delta) = save_calloc("delta", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 847, (delta_xsize), sizeof(*(delta))); |
848 | for (j = 0; j < delta_xsize; j++) |
849 | { |
850 | snew(delta[j], del_lev+1)(delta[j]) = save_calloc("delta[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 850, (del_lev+1), sizeof(*(delta[j]))); |
851 | } |
852 | if (avl > 0) |
853 | { |
854 | snew(rmsdav_mat, tel_mat)(rmsdav_mat) = save_calloc("rmsdav_mat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 854, (tel_mat), sizeof(*(rmsdav_mat))); |
855 | for (j = 0; j < tel_mat; j++) |
856 | { |
857 | snew(rmsdav_mat[j], tel_mat)(rmsdav_mat[j]) = save_calloc("rmsdav_mat[j]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 857, (tel_mat), sizeof(*(rmsdav_mat[j]))); |
858 | } |
859 | } |
860 | } |
861 | |
862 | if (bFitAll) |
863 | { |
864 | snew(mat_x2_j, natoms)(mat_x2_j) = save_calloc("mat_x2_j", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 864, (natoms), sizeof(*(mat_x2_j))); |
865 | } |
866 | for (i = 0; i < tel_mat; i++) |
867 | { |
868 | axis[i] = time[freq*i]; |
869 | fprintf(stderrstderr, "\r element %5d; time %5.2f ", i, axis[i]); |
870 | if (bMat) |
871 | { |
872 | snew(rmsd_mat[i], tel_mat2)(rmsd_mat[i]) = save_calloc("rmsd_mat[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 872, (tel_mat2), sizeof(*(rmsd_mat[i]))); |
873 | } |
874 | if (bBond) |
875 | { |
876 | snew(bond_mat[i], tel_mat2)(bond_mat[i]) = save_calloc("bond_mat[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 876, (tel_mat2), sizeof(*(bond_mat[i]))); |
877 | } |
878 | for (j = 0; j < tel_mat2; j++) |
879 | { |
880 | if (bFitAll) |
881 | { |
882 | for (k = 0; k < n_ind_m; k++) |
883 | { |
884 | copy_rvec(mat_x2[j][k], mat_x2_j[k]); |
885 | } |
886 | do_fit(n_ind_m, w_rls_m, mat_x[i], mat_x2_j); |
887 | } |
888 | else |
889 | { |
890 | mat_x2_j = mat_x2[j]; |
891 | } |
892 | if (bMat) |
893 | { |
894 | if (bFile2 || (i < j)) |
895 | { |
896 | rmsd_mat[i][j] = |
897 | calc_similar_ind(ewhat != ewRMSD, irms[0], ind_rms_m, |
898 | w_rms_m, mat_x[i], mat_x2_j); |
899 | if (rmsd_mat[i][j] > rmsd_max) |
900 | { |
901 | rmsd_max = rmsd_mat[i][j]; |
902 | } |
903 | if (rmsd_mat[i][j] < rmsd_min) |
904 | { |
905 | rmsd_min = rmsd_mat[i][j]; |
906 | } |
907 | rmsd_avg += rmsd_mat[i][j]; |
908 | } |
909 | else |
910 | { |
911 | rmsd_mat[i][j] = rmsd_mat[j][i]; |
912 | } |
913 | } |
914 | if (bBond) |
915 | { |
916 | if (bFile2 || (i <= j)) |
917 | { |
918 | ang = 0.0; |
919 | for (m = 0; m < ibond; m++) |
920 | { |
921 | rvec_sub(mat_x[i][ind_bond1[m]], mat_x[i][ind_bond2[m]], vec1); |
922 | rvec_sub(mat_x2_j[ind_bond1[m]], mat_x2_j[ind_bond2[m]], vec2); |
923 | ang += acos(cos_angle(vec1, vec2)); |
924 | } |
925 | bond_mat[i][j] = ang*180.0/(M_PI3.14159265358979323846*ibond); |
926 | if (bond_mat[i][j] > bond_max) |
927 | { |
928 | bond_max = bond_mat[i][j]; |
929 | } |
930 | if (bond_mat[i][j] < bond_min) |
931 | { |
932 | bond_min = bond_mat[i][j]; |
933 | } |
934 | } |
935 | else |
936 | { |
937 | bond_mat[i][j] = bond_mat[j][i]; |
938 | } |
939 | } |
940 | } |
941 | } |
942 | if (bFile2) |
943 | { |
944 | rmsd_avg /= tel_mat*tel_mat2; |
945 | } |
946 | else |
947 | { |
948 | rmsd_avg /= tel_mat*(tel_mat - 1)/2; |
949 | } |
950 | if (bMat && (avl > 0)) |
951 | { |
952 | rmsd_max = 0.0; |
953 | rmsd_min = 0.0; |
954 | rmsd_avg = 0.0; |
955 | for (j = 0; j < tel_mat-1; j++) |
956 | { |
957 | for (i = j+1; i < tel_mat; i++) |
958 | { |
959 | av_tot = 0; |
960 | weight_tot = 0; |
961 | for (my = -avl; my <= avl; my++) |
962 | { |
963 | if ((j+my >= 0) && (j+my < tel_mat)) |
964 | { |
965 | abs_my = abs(my); |
966 | for (mx = -avl; mx <= avl; mx++) |
967 | { |
968 | if ((i+mx >= 0) && (i+mx < tel_mat)) |
969 | { |
970 | weight = (real)(avl+1-max(abs(mx), abs_my)(((abs(mx)) > (abs_my)) ? (abs(mx)) : (abs_my) )); |
971 | av_tot += weight*rmsd_mat[i+mx][j+my]; |
972 | weight_tot += weight; |
973 | } |
974 | } |
975 | } |
976 | } |
977 | rmsdav_mat[i][j] = av_tot/weight_tot; |
978 | rmsdav_mat[j][i] = rmsdav_mat[i][j]; |
979 | if (rmsdav_mat[i][j] > rmsd_max) |
980 | { |
981 | rmsd_max = rmsdav_mat[i][j]; |
982 | } |
983 | } |
984 | } |
985 | rmsd_mat = rmsdav_mat; |
986 | } |
987 | |
988 | if (bMat) |
989 | { |
990 | fprintf(stderrstderr, "\n%s: Min %f, Max %f, Avg %f\n", |
991 | whatname[ewhat], rmsd_min, rmsd_max, rmsd_avg); |
992 | rlo.r = 1; rlo.g = 1; rlo.b = 1; |
993 | rhi.r = 0; rhi.g = 0; rhi.b = 0; |
994 | if (rmsd_user_max != -1) |
995 | { |
996 | rmsd_max = rmsd_user_max; |
997 | } |
998 | if (rmsd_user_min != -1) |
999 | { |
1000 | rmsd_min = rmsd_user_min; |
1001 | } |
1002 | if ((rmsd_user_max != -1) || (rmsd_user_min != -1)) |
1003 | { |
1004 | fprintf(stderrstderr, "Min and Max value set to resp. %f and %f\n", |
1005 | rmsd_min, rmsd_max); |
1006 | } |
1007 | sprintf(buf, "%s %s matrix", gn_rms[0], whatname[ewhat]); |
1008 | write_xpm(opt2FILE("-m", NFILE, fnm, "w")gmx_ffopen(opt2fn("-m", ((int)(sizeof(fnm)/sizeof((fnm)[0]))) , fnm), "w"), 0, buf, whatlabel[ewhat], |
1009 | output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2, |
1010 | axis, axis2, rmsd_mat, rmsd_min, rmsd_max, rlo, rhi, &nlevels); |
1011 | /* Print the distribution of RMSD values */ |
1012 | if (opt2bSet("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
1013 | { |
1014 | low_rmsd_dist(opt2fn("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), rmsd_max, tel_mat, rmsd_mat, oenv); |
1015 | } |
1016 | |
1017 | if (bDelta) |
1018 | { |
1019 | snew(delta_tot, delta_xsize)(delta_tot) = save_calloc("delta_tot", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 1019, (delta_xsize), sizeof(*(delta_tot))); |
1020 | for (j = 0; j < tel_mat-1; j++) |
1021 | { |
1022 | for (i = j+1; i < tel_mat; i++) |
1023 | { |
1024 | mx = i-j; |
1025 | if (mx < tel_mat/2) |
1026 | { |
1027 | if (bDeltaLog) |
1028 | { |
1029 | mx = (int)(log(mx)*delta_scalex+0.5); |
1030 | } |
1031 | my = (int)(rmsd_mat[i][j]*delta_scaley*del_lev+0.5); |
1032 | delta_tot[mx] += 1.0; |
1033 | if ((rmsd_mat[i][j] >= 0) && (rmsd_mat[i][j] <= delta_maxy)) |
1034 | { |
1035 | delta[mx][my] += 1.0; |
1036 | } |
1037 | } |
1038 | } |
1039 | } |
1040 | delta_max = 0; |
1041 | for (i = 0; i < delta_xsize; i++) |
1042 | { |
1043 | if (delta_tot[i] > 0.0) |
1044 | { |
1045 | delta_tot[i] = 1.0/delta_tot[i]; |
1046 | for (j = 0; j <= del_lev; j++) |
1047 | { |
1048 | delta[i][j] *= delta_tot[i]; |
1049 | if (delta[i][j] > delta_max) |
1050 | { |
1051 | delta_max = delta[i][j]; |
1052 | } |
1053 | } |
1054 | } |
1055 | } |
1056 | fprintf(stderrstderr, "Maximum in delta matrix: %f\n", delta_max); |
1057 | snew(del_xaxis, delta_xsize)(del_xaxis) = save_calloc("del_xaxis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 1057, (delta_xsize), sizeof(*(del_xaxis))); |
1058 | snew(del_yaxis, del_lev+1)(del_yaxis) = save_calloc("del_yaxis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 1058, (del_lev+1), sizeof(*(del_yaxis))); |
1059 | for (i = 0; i < delta_xsize; i++) |
1060 | { |
1061 | del_xaxis[i] = axis[i]-axis[0]; |
1062 | } |
1063 | for (i = 0; i < del_lev+1; i++) |
1064 | { |
1065 | del_yaxis[i] = delta_maxy*i/del_lev; |
1066 | } |
1067 | sprintf(buf, "%s %s vs. delta t", gn_rms[0], whatname[ewhat]); |
1068 | fp = gmx_ffopen("delta.xpm", "w"); |
1069 | write_xpm(fp, 0, buf, "density", output_env_get_time_label(oenv), whatlabel[ewhat], |
1070 | delta_xsize, del_lev+1, del_xaxis, del_yaxis, |
1071 | delta, 0.0, delta_max, rlo, rhi, &nlevels); |
1072 | gmx_ffclose(fp); |
1073 | } |
1074 | if (opt2bSet("-bin", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm)) |
1075 | { |
1076 | /* NB: File must be binary if we use fwrite */ |
1077 | fp = ftp2FILE(efDAT, NFILE, fnm, "wb")gmx_ffopen(ftp2fn(efDAT, ((int)(sizeof(fnm)/sizeof((fnm)[0])) ), fnm), "wb"); |
1078 | for (i = 0; i < tel_mat; i++) |
1079 | { |
1080 | if (fwrite(rmsd_mat[i], sizeof(**rmsd_mat), tel_mat2, fp) != tel_mat2) |
1081 | { |
1082 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_rms.c" , 1082, "Error writing to output file"); |
1083 | } |
1084 | } |
1085 | gmx_ffclose(fp); |
1086 | } |
1087 | } |
1088 | if (bBond) |
1089 | { |
1090 | fprintf(stderrstderr, "\nMin. angle: %f, Max. angle: %f\n", bond_min, bond_max); |
1091 | if (bond_user_max != -1) |
1092 | { |
1093 | bond_max = bond_user_max; |
1094 | } |
1095 | if (bond_user_min != -1) |
1096 | { |
1097 | bond_min = bond_user_min; |
1098 | } |
1099 | if ((bond_user_max != -1) || (bond_user_min != -1)) |
1100 | { |
1101 | fprintf(stderrstderr, "Bond angle Min and Max set to:\n" |
1102 | "Min. angle: %f, Max. angle: %f\n", bond_min, bond_max); |
1103 | } |
1104 | rlo.r = 1; rlo.g = 1; rlo.b = 1; |
1105 | rhi.r = 0; rhi.g = 0; rhi.b = 0; |
1106 | sprintf(buf, "%s av. bond angle deviation", gn_rms[0]); |
1107 | write_xpm(opt2FILE("-bm", NFILE, fnm, "w")gmx_ffopen(opt2fn("-bm", ((int)(sizeof(fnm)/sizeof((fnm)[0])) ), fnm), "w"), 0, buf, "degrees", |
1108 | output_env_get_time_label(oenv), output_env_get_time_label(oenv), tel_mat, tel_mat2, |
1109 | axis, axis2, bond_mat, bond_min, bond_max, rlo, rhi, &nlevels); |
1110 | } |
1111 | } |
1112 | |
1113 | bAv = opt2bSet("-a", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm); |
1114 | |
1115 | /* Write the RMSD's to file */ |
1116 | if (!bPrev) |
1117 | { |
1118 | sprintf(buf, "%s", whatxvgname[ewhat]); |
1119 | } |
1120 | else |
1121 | { |
1122 | sprintf(buf, "%s with frame %g %s ago", whatxvgname[ewhat], |
1123 | time[prev*freq]-time[0], output_env_get_time_label(oenv)); |
1124 | } |
1125 | fp = xvgropen(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), buf, output_env_get_xvgr_tlabel(oenv), |
1126 | whatxvglabel[ewhat], oenv); |
1127 | if (output_env_get_print_xvgr_codes(oenv)) |
1128 | { |
1129 | fprintf(fp, "@ subtitle \"%s%s after %s%s%s\"\n", |
1130 | (nrms == 1) ? "" : "of ", gn_rms[0], fitgraphlabel[efit], |
1131 | bFit ? " to " : "", bFit ? gn_fit : ""); |
1132 | } |
1133 | if (nrms != 1) |
1134 | { |
1135 | xvgr_legend(fp, nrms, (const char**)gn_rms, oenv); |
1136 | } |
1137 | for (i = 0; (i < teller); i++) |
1138 | { |
1139 | if (bSplit && i > 0 && |
1140 | abs(time[bPrev ? freq*i : i]/output_env_get_time_factor(oenv)) < 1e-5) |
1141 | { |
1142 | fprintf(fp, "&\n"); |
1143 | } |
1144 | fprintf(fp, "%12.7f", time[bPrev ? freq*i : i]); |
1145 | for (j = 0; (j < nrms); j++) |
1146 | { |
1147 | fprintf(fp, " %12.7f", rls[j][i]); |
1148 | if (bAv) |
1149 | { |
1150 | rlstot += rls[j][i]; |
1151 | } |
1152 | } |
1153 | fprintf(fp, "\n"); |
1154 | } |
1155 | gmx_ffclose(fp); |
1156 | |
1157 | if (bMirror) |
1158 | { |
1159 | /* Write the mirror RMSD's to file */ |
1160 | sprintf(buf, "%s with Mirror", whatxvgname[ewhat]); |
1161 | sprintf(buf2, "Mirror %s", whatxvglabel[ewhat]); |
1162 | fp = xvgropen(opt2fn("-mir", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), buf, output_env_get_xvgr_tlabel(oenv), |
1163 | buf2, oenv); |
1164 | if (nrms == 1) |
1165 | { |
1166 | if (output_env_get_print_xvgr_codes(oenv)) |
1167 | { |
1168 | fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n", |
1169 | gn_rms[0], gn_fit); |
1170 | } |
1171 | } |
1172 | else |
1173 | { |
1174 | if (output_env_get_print_xvgr_codes(oenv)) |
1175 | { |
1176 | fprintf(fp, "@ subtitle \"after lsq fit to mirror %s\"\n", gn_fit); |
1177 | } |
1178 | xvgr_legend(fp, nrms, (const char**)gn_rms, oenv); |
1179 | } |
1180 | for (i = 0; (i < teller); i++) |
1181 | { |
1182 | if (bSplit && i > 0 && abs(time[i]) < 1e-5) |
1183 | { |
1184 | fprintf(fp, "&\n"); |
1185 | } |
1186 | fprintf(fp, "%12.7f", time[i]); |
1187 | for (j = 0; (j < nrms); j++) |
1188 | { |
1189 | fprintf(fp, " %12.7f", rlsm[j][i]); |
1190 | } |
1191 | fprintf(fp, "\n"); |
1192 | } |
1193 | gmx_ffclose(fp); |
1194 | } |
1195 | |
1196 | if (bAv) |
1197 | { |
1198 | sprintf(buf, "Average %s", whatxvgname[ewhat]); |
1199 | sprintf(buf2, "Average %s", whatxvglabel[ewhat]); |
1200 | fp = xvgropen(opt2fn("-a", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), buf, "Residue", buf2, oenv); |
1201 | for (j = 0; (j < nrms); j++) |
1202 | { |
1203 | fprintf(fp, "%10d %10g\n", j, rlstot/teller); |
1204 | } |
1205 | gmx_ffclose(fp); |
1206 | } |
1207 | |
1208 | if (bNorm) |
1209 | { |
1210 | fp = xvgropen("aver.xvg", gn_rms[0], "Residue", whatxvglabel[ewhat], oenv); |
1211 | for (j = 0; (j < irms[0]); j++) |
1212 | { |
1213 | fprintf(fp, "%10d %10g\n", j, rlsnorm[j]/teller); |
1214 | } |
1215 | gmx_ffclose(fp); |
1216 | } |
1217 | do_view(oenv, opt2fn_null("-a", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "-graphtype bar"); |
1218 | do_view(oenv, opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
1219 | do_view(oenv, opt2fn_null("-mir", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
1220 | do_view(oenv, opt2fn_null("-m", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
1221 | do_view(oenv, opt2fn_null("-bm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
1222 | do_view(oenv, opt2fn_null("-dist", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
1223 | |
1224 | return 0; |
1225 | } |