Bug Summary

File:gromacs/gmxana/gmx_rms.c
Location:line 1168, column 17
Description:Function call argument is an uninitialized value

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 <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
65static 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
98int 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;
1
'gn_fit' declared without an initial value
→
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)
2
←
Taking false branch
→
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)
3
←
Assuming 'ewhat' is not equal to ewRho
→
4
←
Assuming 'ewhat' is not equal to ewRhoSc
→
5
←
Taking false branch
→
282 {
283 please_cite(stdoutstdout, "Maiorov95");
284 }
285 efit = nenum(fit);
286 bFit = efit == efFit;
6
←
Assuming 'efit' is not equal to efFit
→
287 bReset = efit == efReset;
7
←
Assuming 'efit' is not equal to efReset
→
288 if (bFit)
8
←
Taking false branch
→
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)
9
←
Taking false branch
→
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)
10
←
Assuming 'bFreq2' is not equal to 0
→
11
←
Taking false branch
→
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)
12
←
Taking false branch
→
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)
13
←
Assuming 'bDelta' is 0
→
14
←
Taking false branch
→
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)
15
←
Assuming 'bTop' is not equal to 0
→
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)
16
←
Taking false branch
→
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)
17
←
Taking false branch
→
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)
18
←
Taking false branch
→
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]);
19
←
'?' condition is false
→
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)
20
←
Assuming 'bNorm' is 0
→
21
←
Taking false branch
→
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++)
22
←
Loop condition is true. Entering loop body
→
23
←
Loop condition is false. Execution continues on line 435
→
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)
24
←
Assuming 'bMirror' is not equal to 0
→
25
←
Taking true branch
→
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++)
26
←
Loop condition is true. Entering loop body
→
27
←
Loop condition is false. Execution continues on line 443
→
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++)
28
←
Loop condition is true. Entering loop body
→
32
←
Loop condition is false. Execution continues on line 469
→
445 {
446 bMass = FALSE0;
447 for (i = 0; i < irms[j]; i++)
29
←
Loop condition is false. Execution continues on line 459
→
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)
30
←
Taking true branch
→
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++)
31
←
Loop condition is false. Execution continues on line 444
→
463 {
464 w_rms[ind_rms[j][i]] = 1;
465 }
466 }
467 }
468 /* Prepare reference frame */
469 if (bPBC)
33
←
Taking true branch
→
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)
34
←
Taking false branch
→
475 {
476 reset_x(ifit, ind_fit, top.atoms.nr, NULL((void*)0), xp, w_rls);
477 }
478 if (bMirror)
35
←
Taking true branch
→
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++)
36
←
Loop condition is false. Execution continues on line 488
→
483 {
484 copy_rvec(xp[i], xm[i]);
485 xm[i][XX0] = -xm[i][XX0];
486 }
487 }
488 if (ewhat == ewRhoSc)
37
←
Taking false branch
→
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)
38
←
Taking false branch
→
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)
39
←
Taking false branch
→
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)
40
←
Taking false branch
→
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;
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
57
←
Loop condition is false. Exiting loop
→
613 {
614 if (bPBC)
41
←
Taking true branch
→
615 {
616 gmx_rmpbc(gpbc, natoms, box, x);
617 }
618
619 if (bReset)
42
←
Taking false branch
→
620 {
621 reset_x(ifit, ind_fit, natoms, NULL((void*)0), x, w_rls);
622 }
623 if (ewhat == ewRhoSc)
43
←
Taking false branch
→
624 {
625 norm_princ(&top.atoms, ifit, ind_fit, natoms, x);
626 }
627
628 if (bFit)
44
←
Taking false branch
→
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)
45
←
Taking true branch
→
635 {
636 /* keep frame for matrix calculation */
637 if (bMat || bBond || bPrev)
46
←
Taking false branch
→
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)
47
←
Taking false branch
→
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++)
48
←
Loop condition is true. Entering loop body
→
50
←
Loop condition is false. Execution continues on line 678
→
674 {
675 rls[j][teller] =
676 calc_similar_ind(ewhat != ewRMSD, irms[j], ind_rms[j], w_rms, x, xp);
49
←
Assuming 'ewhat' is equal to ewRMSD
→
677 }
678 if (bNorm)
51
←
Taking false branch
→
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)
52
←
Taking true branch
→
688 {
689 if (bFit)
53
←
Taking false branch
→
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++)
54
←
Loop condition is true. Entering loop body
→
55
←
Loop condition is false. Execution continues on line 701
→
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)
56
←
Taking false branch
→
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)
58
←
Taking false branch
→
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)
59
←
Taking false branch
→
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)
60
←
Taking true branch
→
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))
61
←
Taking false branch
→
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)
62
←
Taking false branch
→
1134 {
1135 xvgr_legend(fp, nrms, (const char**)gn_rms, oenv);
1136 }
1137 for (i = 0; (i < teller); i++)
63
←
Loop condition is true. Entering loop body
→
69
←
Loop condition is false. Execution continues on line 1155
→
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]);
64
←
'?' condition is false
→
1145 for (j = 0; (j < nrms); j++)
65
←
Loop condition is true. Entering loop body
→
68
←
Loop condition is false. Execution continues on line 1153
→
1146 {
1147 fprintf(fp, " %12.7f", rls[j][i]);
1148 if (bAv)
66
←
Assuming 'bAv' is 0
→
67
←
Taking false branch
→
1149 {
1150 rlstot += rls[j][i];
1151 }
1152 }
1153 fprintf(fp, "\n");
1154 }
1155 gmx_ffclose(fp);
1156
1157 if (bMirror)
70
←
Taking true branch
→
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)
71
←
Taking true branch
→
1165 {
1166 if (output_env_get_print_xvgr_codes(oenv))
72
←
Taking true branch
→
1167 {
1168 fprintf(fp, "@ subtitle \"of %s after lsq fit to mirror of %s\"\n",
73
←
Function call argument is an uninitialized value
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}