3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
43 #ifdef HAVE_SYS_TIME_H
70 #include "gromacs/linearalgebra/eigensolver.h"
72 /* Portable version of ctime_r implemented in src/gmxlib/string2.c, but we do not want it declared in public installed headers */
74 gmx_ctime_r(const time_t *clock, char *buf, int n);
77 int gmx_covar(int argc, char *argv[])
79 const char *desc[] = {
80 "[TT]g_covar[tt] calculates and diagonalizes the (mass-weighted)",
82 "All structures are fitted to the structure in the structure file.",
83 "When this is not a run input file periodicity will not be taken into",
84 "account. When the fit and analysis groups are identical and the analysis",
85 "is non mass-weighted, the fit will also be non mass-weighted.",
87 "The eigenvectors are written to a trajectory file ([TT]-v[tt]).",
88 "When the same atoms are used for the fit and the covariance analysis,",
89 "the reference structure for the fit is written first with t=-1.",
90 "The average (or reference when [TT]-ref[tt] is used) structure is",
91 "written with t=0, the eigenvectors",
92 "are written as frames with the eigenvector number as timestamp.",
94 "The eigenvectors can be analyzed with [TT]g_anaeig[tt].",
96 "Option [TT]-ascii[tt] writes the whole covariance matrix to",
97 "an ASCII file. The order of the elements is: x1x1, x1y1, x1z1, x1x2, ...",
99 "Option [TT]-xpm[tt] writes the whole covariance matrix to an [TT].xpm[tt] file.",
101 "Option [TT]-xpma[tt] writes the atomic covariance matrix to an [TT].xpm[tt] file,",
102 "i.e. for each atom pair the sum of the xx, yy and zz covariances is",
105 "Note that the diagonalization of a matrix requires memory and time",
106 "that will increase at least as fast as than the square of the number",
107 "of atoms involved. It is easy to run out of memory, in which",
108 "case this tool will probably exit with a 'Segmentation fault'. You",
109 "should consider carefully whether a reduced set of atoms will meet",
110 "your needs for lower costs."
112 static gmx_bool bFit = TRUE, bRef = FALSE, bM = FALSE, bPBC = TRUE;
115 { "-fit", FALSE, etBOOL, {&bFit},
116 "Fit to a reference structure"},
117 { "-ref", FALSE, etBOOL, {&bRef},
118 "Use the deviation from the conformation in the structure file instead of from the average" },
119 { "-mwa", FALSE, etBOOL, {&bM},
120 "Mass-weighted covariance analysis"},
121 { "-last", FALSE, etINT, {&end},
122 "Last eigenvector to write away (-1 is till the last)" },
123 { "-pbc", FALSE, etBOOL, {&bPBC},
124 "Apply corrections for periodic boundary conditions" }
132 rvec *x, *xread, *xref, *xav, *xproj;
134 real *sqrtm, *mat, *eigenvalues, sum, trace, inv_nframes;
135 real t, tstart, tend, **mat2;
136 real xj, *w_rls = NULL;
137 real min, max, *axis;
139 int natoms, nat, count, nframes0, nframes, nlevels;
140 gmx_large_int_t ndim, i, j, k, l;
142 const char *fitfile, *trxfile, *ndxfile;
143 const char *eigvalfile, *eigvecfile, *averfile, *logfile;
144 const char *asciifile, *xpmfile, *xpmafile;
145 char str[STRLEN], *fitname, *ananame, *pcwd;
147 atom_id *index, *ifit;
148 gmx_bool bDiffMass1, bDiffMass2;
150 char timebuf[STRLEN];
154 gmx_rmpbc_t gpbc = NULL;
157 { efTRX, "-f", NULL, ffREAD },
158 { efTPS, NULL, NULL, ffREAD },
159 { efNDX, NULL, NULL, ffOPTRD },
160 { efXVG, NULL, "eigenval", ffWRITE },
161 { efTRN, "-v", "eigenvec", ffWRITE },
162 { efSTO, "-av", "average.pdb", ffWRITE },
163 { efLOG, NULL, "covar", ffWRITE },
164 { efDAT, "-ascii", "covar", ffOPTWR },
165 { efXPM, "-xpm", "covar", ffOPTWR },
166 { efXPM, "-xpma", "covara", ffOPTWR }
168 #define NFILE asize(fnm)
170 parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_BE_NICE,
171 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv);
175 fitfile = ftp2fn(efTPS, NFILE, fnm);
176 trxfile = ftp2fn(efTRX, NFILE, fnm);
177 ndxfile = ftp2fn_null(efNDX, NFILE, fnm);
178 eigvalfile = ftp2fn(efXVG, NFILE, fnm);
179 eigvecfile = ftp2fn(efTRN, NFILE, fnm);
180 averfile = ftp2fn(efSTO, NFILE, fnm);
181 logfile = ftp2fn(efLOG, NFILE, fnm);
182 asciifile = opt2fn_null("-ascii", NFILE, fnm);
183 xpmfile = opt2fn_null("-xpm", NFILE, fnm);
184 xpmafile = opt2fn_null("-xpma", NFILE, fnm);
186 read_tps_conf(fitfile, str, &top, &ePBC, &xref, NULL, box, TRUE);
191 printf("\nChoose a group for the least squares fit\n");
192 get_index(atoms, ndxfile, 1, &nfit, &ifit, &fitname);
195 gmx_fatal(FARGS, "Need >= 3 points to fit!\n");
202 printf("\nChoose a group for the covariance analysis\n");
203 get_index(atoms, ndxfile, 1, &natoms, &index, &ananame);
208 snew(w_rls, atoms->nr);
209 for (i = 0; (i < nfit); i++)
211 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
214 bDiffMass1 = bDiffMass1 || (w_rls[ifit[i]] != w_rls[ifit[i-1]]);
220 for (i = 0; (i < natoms); i++)
224 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
227 bDiffMass2 = bDiffMass2 || (sqrtm[i] != sqrtm[i-1]);
236 if (bFit && bDiffMass1 && !bDiffMass2)
238 bDiffMass1 = natoms != nfit;
240 for (i = 0; (i < natoms) && !bDiffMass1; i++)
242 bDiffMass1 = index[i] != ifit[i];
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++)
252 w_rls[ifit[i]] = 1.0;
257 /* Prepare reference frame */
260 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr, box);
261 gmx_rmpbc(gpbc, atoms->nr, box, xref);
265 reset_x(nfit, ifit, atoms->nr, NULL, xref, w_rls);
271 if (sqrt(GMX_LARGE_INT_MAX) < ndim)
273 gmx_fatal(FARGS, "Number of degrees of freedoms to large for matrix.\n");
275 snew(mat, ndim*ndim);
277 fprintf(stderr, "Calculating the average structure ...\n");
279 nat = read_first_x(oenv, &status, trxfile, &t, &xread, box);
280 if (nat != atoms->nr)
282 fprintf(stderr, "\nWARNING: number of atoms in tpx (%d) and trajectory (%d) do not match\n", natoms, nat);
287 /* calculate x: a fitted struture of the selected atoms */
290 gmx_rmpbc(gpbc, nat, box, xread);
294 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
295 do_fit(nat, w_rls, xref, xread);
297 for (i = 0; i < natoms; i++)
299 rvec_inc(xav[i], xread[index[i]]);
302 while (read_next_x(oenv, status, &t, nat, xread, box));
305 inv_nframes = 1.0/nframes0;
306 for (i = 0; i < natoms; i++)
308 for (d = 0; d < DIM; d++)
310 xav[i][d] *= inv_nframes;
311 xread[index[i]][d] = xav[i][d];
314 write_sto_conf_indexed(opt2fn("-av", NFILE, fnm), "Average structure",
315 atoms, xread, NULL, epbcNONE, zerobox, natoms, index);
318 fprintf(stderr, "Constructing covariance matrix (%dx%d) ...\n", (int)ndim, (int)ndim);
320 nat = read_first_x(oenv, &status, trxfile, &t, &xread, box);
326 /* calculate x: a (fitted) structure of the selected atoms */
329 gmx_rmpbc(gpbc, nat, box, xread);
333 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
334 do_fit(nat, w_rls, xref, xread);
338 for (i = 0; i < natoms; i++)
340 rvec_sub(xread[index[i]], xref[index[i]], x[i]);
345 for (i = 0; i < natoms; i++)
347 rvec_sub(xread[index[i]], xav[i], x[i]);
351 for (j = 0; j < natoms; j++)
353 for (dj = 0; dj < DIM; dj++)
357 for (i = j; i < natoms; i++)
360 for (d = 0; d < DIM; d++)
362 mat[l+d] += x[i][d]*xj;
368 while (read_next_x(oenv, status, &t, nat, xread, box) &&
369 (bRef || nframes < nframes0));
371 gmx_rmpbc_done(gpbc);
373 fprintf(stderr, "Read %d frames\n", nframes);
377 /* copy the reference structure to the ouput array x */
379 for (i = 0; i < natoms; i++)
381 copy_rvec(xref[index[i]], xproj[i]);
389 /* correct the covariance matrix for the mass */
390 inv_nframes = 1.0/nframes;
391 for (j = 0; j < natoms; j++)
393 for (dj = 0; dj < DIM; dj++)
395 for (i = j; i < natoms; i++)
397 k = ndim*(DIM*j+dj)+DIM*i;
398 for (d = 0; d < DIM; d++)
400 mat[k+d] = mat[k+d]*inv_nframes*sqrtm[i]*sqrtm[j];
406 /* symmetrize the matrix */
407 for (j = 0; j < ndim; j++)
409 for (i = j; i < ndim; i++)
411 mat[ndim*i+j] = mat[ndim*j+i];
416 for (i = 0; i < ndim; i++)
418 trace += mat[i*ndim+i];
420 fprintf(stderr, "\nTrace of the covariance matrix: %g (%snm^2)\n",
421 trace, bM ? "u " : "");
425 out = ffopen(asciifile, "w");
426 for (j = 0; j < ndim; j++)
428 for (i = 0; i < ndim; i += 3)
430 fprintf(out, "%g %g %g\n",
431 mat[ndim*j+i], mat[ndim*j+i+1], mat[ndim*j+i+2]);
442 for (j = 0; j < ndim; j++)
444 mat2[j] = &(mat[ndim*j]);
445 for (i = 0; i <= j; i++)
447 if (mat2[j][i] < min)
451 if (mat2[j][j] > max)
458 for (i = 0; i < ndim; i++)
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 = ffopen(xpmfile, "w");
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);
479 snew(mat2, ndim/DIM);
480 for (i = 0; i < ndim/DIM; i++)
482 snew(mat2[i], ndim/DIM);
484 for (j = 0; j < ndim/DIM; j++)
486 for (i = 0; i <= j; i++)
489 for (d = 0; d < DIM; d++)
491 mat2[j][i] += mat[ndim*(DIM*j+d)+DIM*i+d];
493 if (mat2[j][i] < min)
497 if (mat2[j][j] > max)
501 mat2[i][j] = mat2[j][i];
504 snew(axis, ndim/DIM);
505 for (i = 0; i < ndim/DIM; i++)
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 = ffopen(xpmafile, "w");
514 write_xpm3(out, 0, "Covariance", bM ? "u nm^2" : "nm^2",
515 "atom", "atom", ndim/DIM, ndim/DIM, axis, axis,
516 mat2, min, 0.0, max, rlo, rmi, rhi, &nlevels);
519 for (i = 0; i < ndim/DIM; i++)
527 /* call diagonalization routine */
529 snew(eigenvalues, ndim);
530 snew(eigenvectors, ndim*ndim);
532 memcpy(eigenvectors, mat, ndim*ndim*sizeof(real));
533 fprintf(stderr, "\nDiagonalizing ...\n");
535 eigensolver(eigenvectors, ndim, 0, ndim, eigenvalues, mat);
538 /* now write the output */
541 for (i = 0; i < ndim; i++)
543 sum += eigenvalues[i];
545 fprintf(stderr, "\nSum of the eigenvalues: %g (%snm^2)\n",
546 sum, bM ? "u " : "");
547 if (fabs(trace-sum) > 0.01*trace)
549 fprintf(stderr, "\nWARNING: eigenvalue sum deviates from the trace of the covariance matrix\n");
552 fprintf(stderr, "\nWriting eigenvalues to %s\n", eigvalfile);
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++)
560 fprintf (out, "%10d %g\n", (int)i+1, eigenvalues[ndim-1-i]);
566 if (nframes-1 < ndim)
577 /* misuse lambda: 0/1 mass weighted analysis no/yes */
580 WriteXref = eWXR_YES;
581 for (i = 0; i < nfit; i++)
583 copy_rvec(xref[ifit[i]], x[i]);
593 /* misuse lambda: -1 for no fit */
594 WriteXref = eWXR_NOFIT;
597 write_eigenvectors(eigvecfile, natoms, mat, TRUE, 1, end,
598 WriteXref, x, bDiffMass1, xproj, bM, eigenvalues);
600 out = ffopen(logfile, "w");
603 gmx_ctime_r(&now, timebuf, STRLEN);
604 fprintf(out, "Covariance analysis log, written %s\n", timebuf);
606 fprintf(out, "Program: %s\n", argv[0]);
607 gmx_getcwd(str, STRLEN);
609 fprintf(out, "Working directory: %s\n\n", str);
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));
615 fprintf(out, "Read reference structure for fit from %s\n", fitfile);
619 fprintf(out, "Read index groups from %s\n", ndxfile);
623 fprintf(out, "Analysis group is '%s' (%d atoms)\n", ananame, natoms);
626 fprintf(out, "Fit group is '%s' (%d atoms)\n", fitname, nfit);
630 fprintf(out, "No fit was used\n");
632 fprintf(out, "Analysis is %smass weighted\n", bDiffMass2 ? "" : "non-");
635 fprintf(out, "Fit is %smass weighted\n", bDiffMass1 ? "" : "non-");
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",
640 fprintf(out, "Trace of the covariance matrix after diagonalizing: %g\n\n",
643 fprintf(out, "Wrote %d eigenvalues to %s\n", (int)ndim, eigvalfile);
644 if (WriteXref == eWXR_YES)
646 fprintf(out, "Wrote reference structure to %s\n", eigvecfile);
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);
653 fprintf(stderr, "Wrote the log to %s\n", logfile);