Bug Summary

File:gromacs/gmxana/gmx_covar.c
Location:line 239, column 9
Description:Value stored to 'i' is never read

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 <string.h>
43#include <time.h>
44
45#ifdef HAVE_SYS_TIME_H
46#include <sys/time.h>
47#endif
48
49#include "gromacs/commandline/pargs.h"
50#include "typedefs.h"
51#include "gromacs/utility/smalloc.h"
52#include "macros.h"
53#include "gromacs/math/vec.h"
54#include "pbc.h"
55#include "gromacs/utility/futil.h"
56#include "index.h"
57#include "gromacs/fileio/confio.h"
58#include "gromacs/fileio/trnio.h"
59#include "mshift.h"
60#include "gromacs/fileio/xvgr.h"
61#include "rmpbc.h"
62#include "txtdump.h"
63#include "gromacs/fileio/matio.h"
64#include "eigio.h"
65#include "physics.h"
66#include "gmx_ana.h"
67#include "gromacs/utility/cstringutil.h"
68#include "gromacs/fileio/trxio.h"
69
70#include "gromacs/linearalgebra/eigensolver.h"
71#include "gromacs/math/do_fit.h"
72#include "gromacs/utility/fatalerror.h"
73
74int gmx_covar(int argc, char *argv[])
75{
76 const char *desc[] = {
77 "[THISMODULE] calculates and diagonalizes the (mass-weighted)",
78 "covariance matrix.",
79 "All structures are fitted to the structure in the structure file.",
80 "When this is not a run input file periodicity will not be taken into",
81 "account. When the fit and analysis groups are identical and the analysis",
82 "is non mass-weighted, the fit will also be non mass-weighted.",
83 "[PAR]",
84 "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
85 "When the same atoms are used for the fit and the covariance analysis,",
86 "the reference structure for the fit is written first with t=-1.",
87 "The average (or reference when [TT]-ref[tt] is used) structure is",
88 "written with t=0, the eigenvectors",
89 "are written as frames with the eigenvector number as timestamp.",
90 "[PAR]",
91 "The eigenvectors can be analyzed with [gmx-anaeig].",
92 "[PAR]",
93 "Option [TT]-ascii[tt] writes the whole covariance matrix to",
94 "an ASCII file. The order of the elements is: x1x1, x1y1, x1z1, x1x2, ...",
95 "[PAR]",
96 "Option [TT]-xpm[tt] writes the whole covariance matrix to an [TT].xpm[tt] file.",
97 "[PAR]",
98 "Option [TT]-xpma[tt] writes the atomic covariance matrix to an [TT].xpm[tt] file,",
99 "i.e. for each atom pair the sum of the xx, yy and zz covariances is",
100 "written.",
101 "[PAR]",
102 "Note that the diagonalization of a matrix requires memory and time",
103 "that will increase at least as fast as than the square of the number",
104 "of atoms involved. It is easy to run out of memory, in which",
105 "case this tool will probably exit with a 'Segmentation fault'. You",
106 "should consider carefully whether a reduced set of atoms will meet",
107 "your needs for lower costs."
108 };
109 static gmx_bool bFit = TRUE1, bRef = FALSE0, bM = FALSE0, bPBC = TRUE1;
110 static int end = -1;
111 t_pargs pa[] = {
112 { "-fit", FALSE0, etBOOL, {&bFit},
113 "Fit to a reference structure"},
114 { "-ref", FALSE0, etBOOL, {&bRef},
115 "Use the deviation from the conformation in the structure file instead of from the average" },
116 { "-mwa", FALSE0, etBOOL, {&bM},
117 "Mass-weighted covariance analysis"},
118 { "-last", FALSE0, etINT, {&end},
119 "Last eigenvector to write away (-1 is till the last)" },
120 { "-pbc", FALSE0, etBOOL, {&bPBC},
121 "Apply corrections for periodic boundary conditions" }
122 };
123 FILE *out;
124 t_trxstatus *status;
125 t_trxstatus *trjout;
126 t_topology top;
127 int ePBC;
128 t_atoms *atoms;
129 rvec *x, *xread, *xref, *xav, *xproj;
130 matrix box, zerobox;
131 real *sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
132 real t, tstart, tend, **mat2;
133 real xj, *w_rls = NULL((void*)0);
134 real min, max, *axis;
135 int ntopatoms, step;
136 int natoms, nat, count, nframes0, nframes, nlevels;
137 gmx_int64_t ndim, i, j, k, l;
138 int WriteXref;
139 const char *fitfile, *trxfile, *ndxfile;
140 const char *eigvalfile, *eigvecfile, *averfile, *logfile;
141 const char *asciifile, *xpmfile, *xpmafile;
142 char str[STRLEN4096], *fitname, *ananame, *pcwd;
143 int d, dj, nfit;
144 atom_id *index, *ifit;
145 gmx_bool bDiffMass1, bDiffMass2;
146 time_t now;
147 char timebuf[STRLEN4096];
148 t_rgb rlo, rmi, rhi;
149 real *eigenvectors;
150 output_env_t oenv;
151 gmx_rmpbc_t gpbc = NULL((void*)0);
152
153 t_filenm fnm[] = {
154 { efTRX, "-f", NULL((void*)0), ffREAD1<<1 },
155 { efTPS, NULL((void*)0), NULL((void*)0), ffREAD1<<1 },
156 { efNDX, NULL((void*)0), NULL((void*)0), ffOPTRD(1<<1 | 1<<3) },
157 { efXVG, NULL((void*)0), "eigenval", ffWRITE1<<2 },
158 { efTRN, "-v", "eigenvec", ffWRITE1<<2 },
159 { efSTO, "-av", "average.pdb", ffWRITE1<<2 },
160 { efLOG, NULL((void*)0), "covar", ffWRITE1<<2 },
161 { efDAT, "-ascii", "covar", ffOPTWR(1<<2| 1<<3) },
162 { efXPM, "-xpm", "covar", ffOPTWR(1<<2| 1<<3) },
163 { efXPM, "-xpma", "covara", ffOPTWR(1<<2| 1<<3) }
164 };
165#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
166
167 if (!parse_common_args(&argc, argv, PCA_CAN_TIME((1<<6) | (1<<7) | (1<<14)) | PCA_TIME_UNIT(1<<15) | PCA_BE_NICE(1<<13),
168 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))
169 {
170 return 0;
171 }
172
173 clear_mat(zerobox);
174
175 fitfile = ftp2fn(efTPS, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
176 trxfile = ftp2fn(efTRX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
177 ndxfile = ftp2fn_null(efNDX, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
178 eigvalfile = ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
179 eigvecfile = ftp2fn(efTRN, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
180 averfile = ftp2fn(efSTO, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
181 logfile = ftp2fn(efLOG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
182 asciifile = opt2fn_null("-ascii", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
183 xpmfile = opt2fn_null("-xpm", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
184 xpmafile = opt2fn_null("-xpma", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
185
186 read_tps_conf(fitfile, str, &top, &ePBC, &xref, NULL((void*)0), box, TRUE1);
187 atoms = &top.atoms;
188
189 if (bFit)
190 {
191 printf("\nChoose a group for the least squares fit\n");
192 get_index(atoms, ndxfile, 1, &nfit, &ifit, &fitname);
193 if (nfit < 3)
194 {
195 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 195, "Need >= 3 points to fit!\n");
196 }
197 }
198 else
199 {
200 nfit = 0;
201 }
202 printf("\nChoose a group for the covariance analysis\n");
203 get_index(atoms, ndxfile, 1, &natoms, &index, &ananame);
204
205 bDiffMass1 = FALSE0;
206 if (bFit)
207 {
208 snew(w_rls, atoms->nr)(w_rls) = save_calloc("w_rls", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 208, (atoms->nr), sizeof(*(w_rls)));
209 for (i = 0; (i < nfit); i++)
210 {
211 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
212 if (i)
213 {
214 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i-1]]);
215 }
216 }
217 }
218 bDiffMass2 = FALSE0;
219 snew(sqrtm, natoms)(sqrtm) = save_calloc("sqrtm", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 219, (natoms), sizeof(*(sqrtm)));
220 for (i = 0; (i < natoms); i++)
221 {
222 if (bM)
223 {
224 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
225 if (i)
226 {
227 bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i-1]);
228 }
229 }
230 else
231 {
232 sqrtm[i] = 1.0;
233 }
234 }
235
236 if (bFit && bDiffMass1 && !bDiffMass2)
237 {
238 bDiffMass1 = natoms != nfit;
239 i = 0;
Value stored to 'i' is never read
240 for (i = 0; (i < natoms) && !bDiffMass1; i++)
241 {
242 bDiffMass1 = index[i] != ifit[i];
243 }
244 if (!bDiffMass1)
245 {
246 fprintf(stderrstderr, "\n"
247 "Note: the fit and analysis group are identical,\n"
248 " while the fit is mass weighted and the analysis is not.\n"
249 " Making the fit non mass weighted.\n\n");
250 for (i = 0; (i < nfit); i++)
251 {
252 w_rls[ifit[i]] = 1.0;
253 }
254 }
255 }
256
257 /* Prepare reference frame */
258 if (bPBC)
259 {
260 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
261 gmx_rmpbc(gpbc, atoms->nr, box, xref);
262 }
263 if (bFit)
264 {
265 reset_x(nfit, ifit, atoms->nr, NULL((void*)0), xref, w_rls);
266 }
267
268 snew(x, natoms)(x) = save_calloc("x", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 268, (natoms), sizeof(*(x)));
269 snew(xav, natoms)(xav) = save_calloc("xav", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 269, (natoms), sizeof(*(xav)));
270 ndim = natoms*DIM3;
271 if (sqrt(GMX_INT64_MAX(9223372036854775807L)) < ndim)
272 {
273 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 273, "Number of degrees of freedoms to large for matrix.\n");
274 }
275 snew(mat, ndim*ndim)(mat) = save_calloc("mat", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 275, (ndim*ndim), sizeof(*(mat)));
276
277 fprintf(stderrstderr, "Calculating the average structure ...\n");
278 nframes0 = 0;
279 nat = read_first_x(oenv, &status, trxfile, &t, &xread, box);
280 if (nat != atoms->nr)
281 {
282 fprintf(stderrstderr, "\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n", natoms, nat);
283 }
284 do
285 {
286 nframes0++;
287 /* calculate x: a fitted struture of the selected atoms */
288 if (bPBC)
289 {
290 gmx_rmpbc(gpbc, nat, box, xread);
291 }
292 if (bFit)
293 {
294 reset_x(nfit, ifit, nat, NULL((void*)0), xread, w_rls);
295 do_fit(nat, w_rls, xref, xread);
296 }
297 for (i = 0; i < natoms; i++)
298 {
299 rvec_inc(xav[i], xread[index[i]]);
300 }
301 }
302 while (read_next_x(oenv, status, &t, xread, box));
303 close_trj(status);
304
305 inv_nframes = 1.0/nframes0;
306 for (i = 0; i < natoms; i++)
307 {
308 for (d = 0; d < DIM3; d++)
309 {
310 xav[i][d] *= inv_nframes;
311 xread[index[i]][d] = xav[i][d];
312 }
313 }
314 write_sto_conf_indexed(opt2fn("-av", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Average structure",
315 atoms, xread, NULL((void*)0), epbcNONE, zerobox, natoms, index);
316 sfree(xread)save_free("xread", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 316, (xread));
317
318 fprintf(stderrstderr, "Constructing covariance matrix (%dx%d) ...\n", (int)ndim, (int)ndim);
319 nframes = 0;
320 nat = read_first_x(oenv, &status, trxfile, &t, &xread, box);
321 tstart = t;
322 do
323 {
324 nframes++;
325 tend = t;
326 /* calculate x: a (fitted) structure of the selected atoms */
327 if (bPBC)
328 {
329 gmx_rmpbc(gpbc, nat, box, xread);
330 }
331 if (bFit)
332 {
333 reset_x(nfit, ifit, nat, NULL((void*)0), xread, w_rls);
334 do_fit(nat, w_rls, xref, xread);
335 }
336 if (bRef)
337 {
338 for (i = 0; i < natoms; i++)
339 {
340 rvec_sub(xread[index[i]], xref[index[i]], x[i]);
341 }
342 }
343 else
344 {
345 for (i = 0; i < natoms; i++)
346 {
347 rvec_sub(xread[index[i]], xav[i], x[i]);
348 }
349 }
350
351 for (j = 0; j < natoms; j++)
352 {
353 for (dj = 0; dj < DIM3; dj++)
354 {
355 k = ndim*(DIM3*j+dj);
356 xj = x[j][dj];
357 for (i = j; i < natoms; i++)
358 {
359 l = k+DIM3*i;
360 for (d = 0; d < DIM3; d++)
361 {
362 mat[l+d] += x[i][d]*xj;
363 }
364 }
365 }
366 }
367 }
368 while (read_next_x(oenv, status, &t, xread, box) &&
369 (bRef || nframes < nframes0));
370 close_trj(status);
371 gmx_rmpbc_done(gpbc);
372
373 fprintf(stderrstderr, "Read %d frames\n", nframes);
374
375 if (bRef)
376 {
377 /* copy the reference structure to the ouput array x */
378 snew(xproj, natoms)(xproj) = save_calloc("xproj", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 378, (natoms), sizeof(*(xproj)));
379 for (i = 0; i < natoms; i++)
380 {
381 copy_rvec(xref[index[i]], xproj[i]);
382 }
383 }
384 else
385 {
386 xproj = xav;
387 }
388
389 /* correct the covariance matrix for the mass */
390 inv_nframes = 1.0/nframes;
391 for (j = 0; j < natoms; j++)
392 {
393 for (dj = 0; dj < DIM3; dj++)
394 {
395 for (i = j; i < natoms; i++)
396 {
397 k = ndim*(DIM3*j+dj)+DIM3*i;
398 for (d = 0; d < DIM3; d++)
399 {
400 mat[k+d] = mat[k+d]*inv_nframes*sqrtm[i]*sqrtm[j];
401 }
402 }
403 }
404 }
405
406 /* symmetrize the matrix */
407 for (j = 0; j < ndim; j++)
408 {
409 for (i = j; i < ndim; i++)
410 {
411 mat[ndim*i+j] = mat[ndim*j+i];
412 }
413 }
414
415 trace = 0;
416 for (i = 0; i < ndim; i++)
417 {
418 trace += mat[i*ndim+i];
419 }
420 fprintf(stderrstderr, "\nTrace of the covariance matrix: %g (%snm^2)\n",
421 trace, bM ? "u " : "");
422
423 if (asciifile)
424 {
425 out = gmx_ffopen(asciifile, "w");
426 for (j = 0; j < ndim; j++)
427 {
428 for (i = 0; i < ndim; i += 3)
429 {
430 fprintf(out, "%g %g %g\n",
431 mat[ndim*j+i], mat[ndim*j+i+1], mat[ndim*j+i+2]);
432 }
433 }
434 gmx_ffclose(out);
435 }
436
437 if (xpmfile)
438 {
439 min = 0;
440 max = 0;
441 snew(mat2, ndim)(mat2) = save_calloc("mat2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 441, (ndim), sizeof(*(mat2)));
442 for (j = 0; j < ndim; j++)
443 {
444 mat2[j] = &(mat[ndim*j]);
445 for (i = 0; i <= j; i++)
446 {
447 if (mat2[j][i] < min)
448 {
449 min = mat2[j][i];
450 }
451 if (mat2[j][j] > max)
452 {
453 max = mat2[j][i];
454 }
455 }
456 }
457 snew(axis, ndim)(axis) = save_calloc("axis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 457, (ndim), sizeof(*(axis)));
458 for (i = 0; i < ndim; i++)
459 {
460 axis[i] = i+1;
461 }
462 rlo.r = 0; rlo.g = 0; rlo.b = 1;
463 rmi.r = 1; rmi.g = 1; rmi.b = 1;
464 rhi.r = 1; rhi.g = 0; rhi.b = 0;
465 out = gmx_ffopen(xpmfile, "w");
466 nlevels = 80;
467 write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
468 "dim", "dim", ndim, ndim, axis, axis,
469 mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
470 gmx_ffclose(out);
471 sfree(axis)save_free("axis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 471, (axis));
472 sfree(mat2)save_free("mat2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 472, (mat2));
473 }
474
475 if (xpmafile)
476 {
477 min = 0;
478 max = 0;
479 snew(mat2, ndim/DIM)(mat2) = save_calloc("mat2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 479, (ndim/3), sizeof(*(mat2)));
480 for (i = 0; i < ndim/DIM3; i++)
481 {
482 snew(mat2[i], ndim/DIM)(mat2[i]) = save_calloc("mat2[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 482, (ndim/3), sizeof(*(mat2[i])));
483 }
484 for (j = 0; j < ndim/DIM3; j++)
485 {
486 for (i = 0; i <= j; i++)
487 {
488 mat2[j][i] = 0;
489 for (d = 0; d < DIM3; d++)
490 {
491 mat2[j][i] += mat[ndim*(DIM3*j+d)+DIM3*i+d];
492 }
493 if (mat2[j][i] < min)
494 {
495 min = mat2[j][i];
496 }
497 if (mat2[j][j] > max)
498 {
499 max = mat2[j][i];
500 }
501 mat2[i][j] = mat2[j][i];
502 }
503 }
504 snew(axis, ndim/DIM)(axis) = save_calloc("axis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 504, (ndim/3), sizeof(*(axis)));
505 for (i = 0; i < ndim/DIM3; i++)
506 {
507 axis[i] = i+1;
508 }
509 rlo.r = 0; rlo.g = 0; rlo.b = 1;
510 rmi.r = 1; rmi.g = 1; rmi.b = 1;
511 rhi.r = 1; rhi.g = 0; rhi.b = 0;
512 out = gmx_ffopen(xpmafile, "w");
513 nlevels = 80;
514 write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
515 "atom", "atom", ndim/DIM3, ndim/DIM3, axis, axis,
516 mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
517 gmx_ffclose(out);
518 sfree(axis)save_free("axis", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 518, (axis));
519 for (i = 0; i < ndim/DIM3; i++)
520 {
521 sfree(mat2[i])save_free("mat2[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 521, (mat2[i]));
522 }
523 sfree(mat2)save_free("mat2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 523, (mat2));
524 }
525
526
527 /* call diagonalization routine */
528
529 snew(eigenvalues, ndim)(eigenvalues) = save_calloc("eigenvalues", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 529, (ndim), sizeof(*(eigenvalues)));
530 snew(eigenvectors, ndim*ndim)(eigenvectors) = save_calloc("eigenvectors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 530, (ndim*ndim), sizeof(*(eigenvectors)));
531
532 memcpy(eigenvectors, mat, ndim*ndim*sizeof(real));
533 fprintf(stderrstderr, "\nDiagonalizing ...\n");
534 fflush(stderrstderr);
535 eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
536 sfree(eigenvectors)save_free("eigenvectors", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_covar.c"
, 536, (eigenvectors));
537
538 /* now write the output */
539
540 sum = 0;
541 for (i = 0; i < ndim; i++)
542 {
543 sum += eigenvalues[i];
544 }
545 fprintf(stderrstderr, "\nSum of the eigenvalues: %g (%snm^2)\n",
546 sum, bM ? "u " : "");
547 if (fabs(trace-sum) > 0.01*trace)
548 {
549 fprintf(stderrstderr, "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
550 }
551
552 fprintf(stderrstderr, "\nWriting eigenvalues to %s\n", eigvalfile);
553
554 sprintf(str, "(%snm\\S2\\N)", bM ? "u " : "");
555 out = xvgropen(eigvalfile,
556 "Eigenvalues of the covariance matrix",
557 "Eigenvector index", str, oenv);
558 for (i = 0; (i < ndim); i++)
559 {
560 fprintf (out, "%10d %g\n", (int)i+1, eigenvalues[ndim-1-i]);
561 }
562 gmx_ffclose(out);
563
564 if (end == -1)
565 {
566 if (nframes-1 < ndim)
567 {
568 end = nframes-1;
569 }
570 else
571 {
572 end = ndim;
573 }
574 }
575 if (bFit)
576 {
577 /* misuse lambda: 0/1 mass weighted analysis no/yes */
578 if (nfit == natoms)
579 {
580 WriteXref = eWXR_YES;
581 for (i = 0; i < nfit; i++)
582 {
583 copy_rvec(xref[ifit[i]], x[i]);
584 }
585 }
586 else
587 {
588 WriteXref = eWXR_NO;
589 }
590 }
591 else
592 {
593 /* misuse lambda: -1 for no fit */
594 WriteXref = eWXR_NOFIT;
595 }
596
597 write_eigenvectors(eigvecfile, natoms, mat, TRUE1, 1, end,
598 WriteXref, x, bDiffMass1, xproj, bM, eigenvalues);
599
600 out = gmx_ffopen(logfile, "w");
601
602 time(&now);
603 gmx_ctime_r(&now, timebuf, STRLEN4096);
604 fprintf(out, "Covariance analysis log, written %s\n", timebuf);
605
606 fprintf(out, "Program: %s\n", argv[0]);
607 gmx_getcwd(str, STRLEN4096);
608
609 fprintf(out, "Working directory: %s\n\n", str);
610
611 fprintf(out, "Read %d frames from %s (time %g to %g %s)\n", nframes, trxfile,
612 output_env_conv_time(oenv, tstart), output_env_conv_time(oenv, tend), output_env_get_time_unit(oenv));
613 if (bFit)
614 {
615 fprintf(out, "Read reference structure for fit from %s\n", fitfile);
616 }
617 if (ndxfile)
618 {
619 fprintf(out, "Read index groups from %s\n", ndxfile);
620 }
621 fprintf(out, "\n");
622
623 fprintf(out, "Analysis group is '%s' (%d atoms)\n", ananame, natoms);
624 if (bFit)
625 {
626 fprintf(out, "Fit group is '%s' (%d atoms)\n", fitname, nfit);
627 }
628 else
629 {
630 fprintf(out, "No fit was used\n");
631 }
632 fprintf(out, "Analysis is %smass weighted\n", bDiffMass2 ? "" : "non-");
633 if (bFit)
634 {
635 fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
636 }
637 fprintf(out, "Diagonalized the %dx%d covariance matrix\n", (int)ndim, (int)ndim);
638 fprintf(out, "Trace of the covariance matrix before diagonalizing: %g\n",
639 trace);
640 fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n",
641 sum);
642
643 fprintf(out, "Wrote %d eigenvalues to %s\n", (int)ndim, eigvalfile);
644 if (WriteXref == eWXR_YES)
645 {
646 fprintf(out, "Wrote reference structure to %s\n", eigvecfile);
647 }
648 fprintf(out, "Wrote average structure to %s and %s\n", averfile, eigvecfile);
649 fprintf(out, "Wrote eigenvectors %d to %d to %s\n", 1, end, eigvecfile);
650
651 gmx_ffclose(out);
652
653 fprintf(stderrstderr, "Wrote the log to %s\n", logfile);
654
655 return 0;
656}