2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/legacyheaders/typedefs.h"
45 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/legacyheaders/macros.h"
47 #include "gromacs/utility/fatalerror.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/utility/futil.h"
50 #include "gromacs/topology/index.h"
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/fileio/confio.h"
53 #include "gromacs/fileio/tpxio.h"
54 #include "gromacs/fileio/trxio.h"
55 #include "gromacs/fileio/matio.h"
56 #include "gromacs/fileio/xvgr.h"
57 #include "gromacs/legacyheaders/viewit.h"
58 #include "gromacs/pbcutil/rmpbc.h"
59 #include "gromacs/legacyheaders/txtdump.h"
61 #include "gromacs/math/units.h"
64 #include "gromacs/math/do_fit.h"
66 static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
69 double hwkT, w, dS, S = 0;
72 hbar = PLANCK1/(2*M_PI);
73 for (i = 0; (i < n-nskip); i++)
77 lambda = eigval[i]*AMU;
78 w = sqrt(BOLTZMANN*temp/lambda)/NANO;
79 hwkT = (hbar*w)/(BOLTZMANN*temp);
80 dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
84 fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
85 i, w, lambda, hwkT, dS);
90 fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
94 fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
98 static void calc_entropy_schlitter(FILE *fp, int n, int nskip,
99 real *eigval, real temp)
105 double hbar, kt, kteh, S;
107 hbar = PLANCK1/(2*M_PI);
109 kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
112 fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
116 for (i = 0; (i < n-nskip); i++)
118 dd = 1+kteh*eigval[i];
123 fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S);
126 const char *proj_unit;
128 static real tick_spacing(real range, int minticks)
137 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
138 while (range/sp < minticks-1)
146 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
147 const char *title, const char *subtitle,
148 const char *xlabel, const char **ylabel,
149 int n, real *x, real **y, real ***sy,
150 real scale_x, gmx_bool bZero, gmx_bool bSplit,
151 const output_env_t oenv)
155 real min, max, xsp, ysp;
157 out = gmx_ffopen(file, "w");
158 if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
160 fprintf(out, "@ autoscale onread none\n");
162 for (g = 0; g < ngraphs; g++)
168 for (i = 0; i < n; i++)
184 for (s = 0; s < nsetspergraph; s++)
186 for (i = 0; i < n; i++)
188 if (sy[g][s][i] < min)
192 if (sy[g][s][i] > max)
205 min = min-0.1*(max-min);
207 max = max+0.1*(max-min);
208 xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
209 ysp = tick_spacing(max-min, 3);
210 if (output_env_get_print_xvgr_codes(oenv))
212 fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
215 fprintf(out, "@ title \"%s\"\n", title);
218 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
223 fprintf(out, "@ xaxis label \"%s\"\n", xlabel);
227 fprintf(out, "@ xaxis ticklabel off\n");
231 fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
232 fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
233 fprintf(out, "@ world ymin %g\n", min);
234 fprintf(out, "@ world ymax %g\n", max);
236 fprintf(out, "@ view xmin 0.15\n");
237 fprintf(out, "@ view xmax 0.85\n");
238 fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs);
239 fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs);
240 fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
241 fprintf(out, "@ xaxis tick major %g\n", xsp);
242 fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
243 fprintf(out, "@ xaxis ticklabel start type spec\n");
244 fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp);
245 fprintf(out, "@ yaxis tick major %g\n", ysp);
246 fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
247 fprintf(out, "@ yaxis ticklabel start type spec\n");
248 fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp);
249 if ((min < 0) && (max > 0))
251 fprintf(out, "@ zeroxaxis bar on\n");
252 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
255 for (s = 0; s < nsetspergraph; s++)
257 for (i = 0; i < n; i++)
259 if (bSplit && i > 0 && fabs(x[i]) < 1e-5)
261 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
263 fprintf(out, "%10.4f %10.5f\n",
264 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
266 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
273 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
274 real *eigval1, int neig1, real *eigval2, int neig2)
278 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
282 n = min(n, min(neig1, neig2));
283 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
286 for (i = 0; i < n; i++)
293 eigval1[i] = sqrt(eigval1[i]);
296 for (i = n; i < neig1; i++)
298 trace1 += eigval1[i];
301 for (i = 0; i < n; i++)
308 eigval2[i] = sqrt(eigval2[i]);
311 for (i = n; i < neig2; i++)
313 trace2 += eigval2[i];
316 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
317 if (neig1 != n || neig2 != n)
319 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
320 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
322 fprintf(stdout, "Square root of the traces: %g and %g\n",
323 sqrt(sum1), sqrt(sum2));
326 for (i = 0; i < n; i++)
329 for (j = 0; j < n; j++)
332 for (k = 0; k < natoms; k++)
334 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
336 tmp += eigval2[j]*ip*ip;
338 sab += eigval1[i]*tmp;
341 samsb2 = sum1+sum2-2*sab;
347 fprintf(stdout, "The overlap of the covariance matrices:\n");
348 fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
349 tmp = 1-sab/sqrt(sum1*sum2);
354 fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
358 static void inprod_matrix(const char *matfile, int natoms,
359 int nvec1, int *eignr1, rvec **eigvec1,
360 int nvec2, int *eignr2, rvec **eigvec2,
361 gmx_bool bSelect, int noutvec, int *outvec)
365 int i, x1, y1, x, y, nlevels;
367 real inp, *t_x, *t_y, max;
375 for (y1 = 0; y1 < nx; y1++)
377 if (outvec[y1] < nvec2)
379 t_y[ny] = eignr2[outvec[y1]]+1;
388 for (y = 0; y < ny; y++)
390 t_y[y] = eignr2[y]+1;
394 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
400 for (x1 = 0; x1 < nx; x1++)
411 t_x[x1] = eignr1[x]+1;
412 fprintf(stderr, " %d", eignr1[x]+1);
413 for (y1 = 0; y1 < ny; y1++)
418 while (outvec[y1] >= nvec2)
428 for (i = 0; i < natoms; i++)
430 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
432 mat[x1][y1] = fabs(inp);
433 if (mat[x1][y1] > max)
439 fprintf(stderr, "\n");
440 rlo.r = 1; rlo.g = 1; rlo.b = 1;
441 rhi.r = 0; rhi.g = 0; rhi.b = 0;
443 out = gmx_ffopen(matfile, "w");
444 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
445 nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
449 static void overlap(const char *outfile, int natoms,
451 int nvec2, int *eignr2, rvec **eigvec2,
452 int noutvec, int *outvec,
453 const output_env_t oenv)
459 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
460 for (i = 0; i < noutvec; i++)
462 fprintf(stderr, "%d ", outvec[i]+1);
464 fprintf(stderr, "\n");
466 out = xvgropen(outfile, "Subspace overlap",
467 "Eigenvectors of trajectory 2", "Overlap", oenv);
468 if (output_env_get_print_xvgr_codes(oenv))
470 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
473 for (x = 0; x < nvec2; x++)
475 for (v = 0; v < noutvec; v++)
479 for (i = 0; i < natoms; i++)
481 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
485 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
491 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
492 const char *projfile, const char *twodplotfile,
493 const char *threedplotfile, const char *filterfile, int skip,
494 const char *extremefile, gmx_bool bExtrAll, real extreme,
495 int nextr, t_atoms *atoms, int natoms, atom_id *index,
496 gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
497 real *sqrtm, rvec *xav,
498 int *eignr, rvec **eigvec,
499 int noutvec, int *outvec, gmx_bool bSplit,
500 const output_env_t oenv)
502 FILE *xvgrout = NULL;
503 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
504 t_trxstatus *out = NULL;
506 int noutvec_extr, imin, imax;
511 real t, inp, **inprod = NULL, min = 0, max = 0;
512 char str[STRLEN], str2[STRLEN], **ylabel, *c;
514 gmx_rmpbc_t gpbc = NULL;
520 noutvec_extr = noutvec;
530 snew(inprod, noutvec+1);
534 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
536 for (i = 0; i < noutvec; i++)
538 fprintf(stderr, "%d ", outvec[i]+1);
540 fprintf(stderr, "\n");
541 out = open_trx(filterfile, "w");
546 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
549 gmx_fatal(FARGS, "the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)", nat, atoms->nr);
555 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
558 for (i = 0; i < nat; i++)
568 gmx_rmpbc(gpbc, nat, box, xread);
570 if (nframes >= snew_size)
573 for (i = 0; i < noutvec+1; i++)
575 srenew(inprod[i], snew_size);
578 inprod[noutvec][nframes] = t;
579 /* calculate x: a fitted struture of the selected atoms */
582 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
583 do_fit(nat, w_rls, xref, xread);
585 for (i = 0; i < natoms; i++)
587 copy_rvec(xread[index[i]], x[i]);
590 for (v = 0; v < noutvec; v++)
593 /* calculate (mass-weighted) projection */
595 for (i = 0; i < natoms; i++)
597 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
598 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
599 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
601 inprod[v][nframes] = inp;
605 for (i = 0; i < natoms; i++)
607 for (d = 0; d < DIM; d++)
609 /* misuse xread for output */
610 xread[index[i]][d] = xav[i][d];
611 for (v = 0; v < noutvec; v++)
613 xread[index[i]][d] +=
614 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
618 write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
624 while (read_next_x(oenv, status, &t, xread, box));
634 snew(xread, atoms->nr);
639 gmx_rmpbc_done(gpbc);
645 snew(ylabel, noutvec);
646 for (v = 0; v < noutvec; v++)
648 sprintf(str, "vec %d", eignr[outvec[v]]+1);
649 ylabel[v] = gmx_strdup(str);
651 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
652 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
653 (const char **)ylabel,
654 nframes, inprod[noutvec], inprod, NULL,
655 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
660 sprintf(str, "projection on eigenvector %d (%s)",
661 eignr[outvec[0]]+1, proj_unit);
662 sprintf(str2, "projection on eigenvector %d (%s)",
663 eignr[outvec[noutvec-1]]+1, proj_unit);
664 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
666 for (i = 0; i < nframes; i++)
668 if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
670 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
672 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
674 gmx_ffclose(xvgrout);
689 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
693 bPDB = fn2ftp(threedplotfile) == efPDB;
695 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
697 b4D = bPDB && (noutvec >= 4);
700 fprintf(stderr, "You have selected four or more eigenvectors:\n"
701 "fourth eigenvector will be plotted "
702 "in bfactor field of pdb file\n");
703 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
704 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
705 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
709 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
710 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
712 init_t_atoms(&atoms, nframes, FALSE);
715 atnm = gmx_strdup("C");
716 resnm = gmx_strdup("PRJ");
720 fact = 10000.0/nframes;
727 for (i = 0; i < nframes; i++)
729 atoms.atomname[i] = &atnm;
730 atoms.atom[i].resind = i;
731 atoms.resinfo[i].name = &resnm;
732 atoms.resinfo[i].nr = ceil(i*fact);
733 atoms.resinfo[i].ic = ' ';
734 x[i][XX] = inprod[0][i];
735 x[i][YY] = inprod[1][i];
736 x[i][ZZ] = inprod[2][i];
742 if ( ( b4D || bSplit ) && bPDB)
744 out = gmx_ffopen(threedplotfile, "w");
745 fprintf(out, "HEADER %s\n", str);
748 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
751 for (i = 0; i < atoms.nr; i++)
753 if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
755 fprintf(out, "TER\n");
758 gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
759 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
762 fprintf(out, "CONECT%5d%5d\n", i, i+1);
766 fprintf(out, "TER\n");
771 write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
773 free_t_atoms(&atoms, FALSE);
778 snew(pmin, noutvec_extr);
779 snew(pmax, noutvec_extr);
782 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
784 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
787 for (v = 0; v < noutvec_extr; v++)
789 for (i = 0; i < nframes; i++)
791 if (inprod[v][i] < inprod[v][imin])
795 if (inprod[v][i] > inprod[v][imax])
800 pmin[v] = inprod[v][imin];
801 pmax[v] = inprod[v][imax];
802 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
804 pmin[v], imin, pmax[v], imax);
812 /* build format string for filename: */
813 strcpy(str, extremefile); /* copy filename */
814 c = strrchr(str, '.'); /* find where extention begins */
815 strcpy(str2, c); /* get extention */
816 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
817 for (v = 0; v < noutvec_extr; v++)
819 /* make filename using format string */
820 if (noutvec_extr == 1)
822 strcpy(str2, extremefile);
826 sprintf(str2, str, eignr[outvec[v]]+1);
828 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
829 nextr, outvec[v]+1, str2);
830 out = open_trx(str2, "w");
831 for (frame = 0; frame < nextr; frame++)
833 if ((extreme == 0) && (nextr <= 3))
835 for (i = 0; i < natoms; i++)
837 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
840 for (i = 0; i < natoms; i++)
842 for (d = 0; d < DIM; d++)
845 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
846 *eigvec[outvec[v]][i][d]/sqrtm[i]);
849 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
856 fprintf(stderr, "\n");
859 static void components(const char *outfile, int natoms,
860 int *eignr, rvec **eigvec,
861 int noutvec, int *outvec,
862 const output_env_t oenv)
866 char str[STRLEN], **ylabel;
868 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
870 snew(ylabel, noutvec);
873 for (i = 0; i < natoms; i++)
877 for (g = 0; g < noutvec; g++)
880 sprintf(str, "vec %d", eignr[v]+1);
881 ylabel[g] = gmx_strdup(str);
883 for (s = 0; s < 4; s++)
885 snew(y[g][s], natoms);
887 for (i = 0; i < natoms; i++)
889 y[g][0][i] = norm(eigvec[v][i]);
890 for (s = 0; s < 3; s++)
892 y[g][s+1][i] = eigvec[v][i][s];
896 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
897 "black: total, red: x, green: y, blue: z",
898 "Atom number", (const char **)ylabel,
899 natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
900 fprintf(stderr, "\n");
903 static void rmsf(const char *outfile, int natoms, real *sqrtm,
904 int *eignr, rvec **eigvec,
905 int noutvec, int *outvec,
906 real *eigval, int neig,
907 const output_env_t oenv)
911 char str[STRLEN], **ylabel;
913 for (i = 0; i < neig; i++)
921 fprintf(stderr, "Writing rmsf to %s\n", outfile);
923 snew(ylabel, noutvec);
926 for (i = 0; i < natoms; i++)
930 for (g = 0; g < noutvec; g++)
933 if (eignr[v] >= neig)
935 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
937 sprintf(str, "vec %d", eignr[v]+1);
938 ylabel[g] = gmx_strdup(str);
940 for (i = 0; i < natoms; i++)
942 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
945 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
946 "Atom number", (const char **)ylabel,
947 natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
948 fprintf(stderr, "\n");
951 int gmx_anaeig(int argc, char *argv[])
953 static const char *desc[] = {
954 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
955 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
956 "([gmx-nmeig]).[PAR]",
958 "When a trajectory is projected on eigenvectors, all structures are",
959 "fitted to the structure in the eigenvector file, if present, otherwise",
960 "to the structure in the structure file. When no run input file is",
961 "supplied, periodicity will not be taken into account. Most analyses",
962 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
963 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
965 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
966 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
968 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
969 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
971 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
972 "[TT]-first[tt] to [TT]-last[tt].",
973 "The projections of a trajectory on the eigenvectors of its",
974 "covariance matrix are called principal components (pc's).",
975 "It is often useful to check the cosine content of the pc's,",
976 "since the pc's of random diffusion are cosines with the number",
977 "of periods equal to half the pc index.",
978 "The cosine content of the pc's can be calculated with the program",
979 "[gmx-analyze].[PAR]",
981 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
982 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
984 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
985 "three selected eigenvectors.[PAR]",
987 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
988 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
990 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
991 "on the average structure and interpolate [TT]-nframes[tt] frames",
992 "between them, or set your own extremes with [TT]-max[tt]. The",
993 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
994 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
995 "will be written to separate files. Chain identifiers will be added",
996 "when writing a [TT].pdb[tt] file with two or three structures (you",
997 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
999 " Overlap calculations between covariance analysis:[BR]",
1000 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1002 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1003 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1004 "in file [TT]-v[tt].[PAR]",
1006 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1007 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1008 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1009 "have been set explicitly.[PAR]",
1011 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1012 "a single number for the overlap between the covariance matrices is",
1013 "generated. The formulas are:[BR]",
1014 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1015 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1016 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1017 "where M1 and M2 are the two covariance matrices and tr is the trace",
1018 "of a matrix. The numbers are proportional to the overlap of the square",
1019 "root of the fluctuations. The normalized overlap is the most useful",
1020 "number, it is 1 for identical matrices and 0 when the sampled",
1021 "subspaces are orthogonal.[PAR]",
1022 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1023 "computed based on the Quasiharmonic approach and based on",
1024 "Schlitter's formula."
1026 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1027 static real max = 0.0, temp = 298.15;
1028 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1030 { "-first", FALSE, etINT, {&first},
1031 "First eigenvector for analysis (-1 is select)" },
1032 { "-last", FALSE, etINT, {&last},
1033 "Last eigenvector for analysis (-1 is till the last)" },
1034 { "-skip", FALSE, etINT, {&skip},
1035 "Only analyse every nr-th frame" },
1036 { "-max", FALSE, etREAL, {&max},
1037 "Maximum for projection of the eigenvector on the average structure, "
1038 "max=0 gives the extremes" },
1039 { "-nframes", FALSE, etINT, {&nextr},
1040 "Number of frames for the extremes output" },
1041 { "-split", FALSE, etBOOL, {&bSplit},
1042 "Split eigenvector projections where time is zero" },
1043 { "-entropy", FALSE, etBOOL, {&bEntropy},
1044 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1045 { "-temp", FALSE, etREAL, {&temp},
1046 "Temperature for entropy calculations" },
1047 { "-nevskip", FALSE, etINT, {&nskip},
1048 "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic approximation. When you do a rotational and/or translational fit prior to the covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which should not be taken into account when computing the entropy." }
1050 #define NPA asize(pa)
1056 t_atoms *atoms = NULL;
1057 rvec *xtop, *xref1, *xref2, *xrefp = NULL;
1058 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1059 int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1060 rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1062 real xid, totmass, *sqrtm, *w_rls, t, lambda;
1065 const char *indexfile;
1068 int nout, *iout, noutvec, *outvec, nfit;
1069 atom_id *index, *ifit;
1070 const char *VecFile, *Vec2File, *topfile;
1071 const char *EigFile, *Eig2File;
1072 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1073 const char *TwoDPlotFile, *ThreeDPlotFile;
1074 const char *FilterFile, *ExtremeFile;
1075 const char *OverlapFile, *InpMatFile;
1076 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1077 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1078 real *eigval1 = NULL, *eigval2 = NULL;
1085 { efTRN, "-v", "eigenvec", ffREAD },
1086 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1087 { efTRX, "-f", NULL, ffOPTRD },
1088 { efTPS, NULL, NULL, ffOPTRD },
1089 { efNDX, NULL, NULL, ffOPTRD },
1090 { efXVG, "-eig", "eigenval", ffOPTRD },
1091 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1092 { efXVG, "-comp", "eigcomp", ffOPTWR },
1093 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1094 { efXVG, "-proj", "proj", ffOPTWR },
1095 { efXVG, "-2d", "2dproj", ffOPTWR },
1096 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1097 { efTRX, "-filt", "filtered", ffOPTWR },
1098 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1099 { efXVG, "-over", "overlap", ffOPTWR },
1100 { efXPM, "-inpr", "inprod", ffOPTWR }
1102 #define NFILE asize(fnm)
1104 if (!parse_common_args(&argc, argv,
1105 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
1106 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1111 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1113 VecFile = opt2fn("-v", NFILE, fnm);
1114 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1115 topfile = ftp2fn(efTPS, NFILE, fnm);
1116 EigFile = opt2fn_null("-eig", NFILE, fnm);
1117 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1118 CompFile = opt2fn_null("-comp", NFILE, fnm);
1119 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1120 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1121 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1122 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1123 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1124 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1125 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1126 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1128 bTop = fn2bTPX(topfile);
1129 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1130 || FilterFile || ExtremeFile;
1132 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1133 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1134 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1135 bVec2 = Vec2File || OverlapFile || InpMatFile;
1136 bM = RmsfFile || bProj;
1137 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1138 || TwoDPlotFile || ThreeDPlotFile;
1139 bIndex = bM || bProj;
1140 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1141 FilterFile || (bIndex && indexfile);
1142 bCompare = Vec2File || Eig2File;
1143 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1145 read_eigenvectors(VecFile, &natoms, &bFit1,
1146 &xref1, &bDMR1, &xav1, &bDMA1,
1147 &nvec1, &eignr1, &eigvec1, &eigval1);
1150 /* Overwrite eigenvalues from separate files if the user provides them */
1151 if (EigFile != NULL)
1153 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1154 if (neig_tmp != neig1)
1156 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1159 srenew(eigval1, neig1);
1160 for (j = 0; j < neig1; j++)
1162 real tmp = eigval1[j];
1163 eigval1[j] = xvgdata[1][j];
1164 if (debug && (eigval1[j] != tmp))
1166 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1167 j, tmp, eigval1[j]);
1170 for (j = 0; j < i; j++)
1175 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1182 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1184 calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1185 calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1192 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1194 read_eigenvectors(Vec2File, &neig2, &bFit2,
1195 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1200 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1204 if (Eig2File != NULL)
1206 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1207 srenew(eigval2, neig2);
1208 for (j = 0; j < neig2; j++)
1210 eigval2[j] = xvgdata[1][j];
1212 for (j = 0; j < i; j++)
1217 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1221 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1225 if ((xref1 == NULL) && (bM || bTraj))
1241 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1242 title, &top, &ePBC, &xtop, NULL, topbox, bM);
1244 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1245 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1246 /* Fitting is only required for the projection */
1251 printf("\nNote: the structure in %s should be the same\n"
1252 " as the one used for the fit in g_covar\n", topfile);
1254 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1255 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1257 snew(w_rls, atoms->nr);
1258 for (i = 0; (i < nfit); i++)
1262 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1266 w_rls[ifit[i]] = 1.0;
1270 snew(xrefp, atoms->nr);
1273 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1276 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d, your selection does not fit the reference structure in the eigenvector file.", nfit, natoms);
1278 for (i = 0; (i < nfit); i++)
1280 copy_rvec(xref1[i], xrefp[ifit[i]]);
1285 /* The top coordinates are the fitting reference */
1286 for (i = 0; (i < nfit); i++)
1288 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1290 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1293 gmx_rmpbc_done(gpbc);
1298 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1299 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1302 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1307 snew(sqrtm, natoms);
1310 proj_unit = "u\\S1/2\\Nnm";
1311 for (i = 0; (i < natoms); i++)
1313 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1319 for (i = 0; (i < natoms); i++)
1329 for (i = 0; (i < natoms); i++)
1331 for (d = 0; (d < DIM); d++)
1333 t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1334 totmass += sqr(sqrtm[i]);
1337 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1338 " %.3f (nm)\n\n", sqrt(t/totmass));
1349 /* make an index from first to last */
1350 nout = last-first+1;
1352 for (i = 0; i < nout; i++)
1354 iout[i] = first-1+i;
1357 else if (ThreeDPlotFile)
1359 /* make an index of first+(0,1,2) and last */
1360 nout = bPDB3D ? 4 : 3;
1361 nout = min(last-first+1, nout);
1369 iout[nout-1] = last-1;
1373 /* make an index of first and last */
1382 printf("Select eigenvectors for output, end your selection with 0\n");
1389 srenew(iout, nout+1);
1390 if (1 != scanf("%d", &iout[nout]))
1392 gmx_fatal(FARGS, "Error reading user input");
1396 while (iout[nout] >= 0);
1400 /* make an index of the eigenvectors which are present */
1403 for (i = 0; i < nout; i++)
1406 while ((j < nvec1) && (eignr1[j] != iout[i]))
1410 if ((j < nvec1) && (eignr1[j] == iout[i]))
1412 outvec[noutvec] = j;
1416 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1419 fprintf(stderr, ":");
1420 for (j = 0; j < noutvec; j++)
1422 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1425 fprintf(stderr, "\n");
1429 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1434 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1440 project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1441 bTop ? &top : NULL, ePBC, topbox,
1442 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1443 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1444 bFit1, xrefp, nfit, ifit, w_rls,
1445 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1451 overlap(OverlapFile, natoms,
1452 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1457 inprod_matrix(InpMatFile, natoms,
1458 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1459 bFirstLastSet, noutvec, outvec);
1464 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1468 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1469 !bCompare && !bEntropy)
1471 fprintf(stderr, "\nIf you want some output,"
1472 " set one (or two or ...) of the output file options\n");
1476 view_all(oenv, NFILE, fnm);