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
45 #include "gmx_fatal.h"
64 static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
67 double hwkT, w, dS, S = 0;
70 hbar = PLANCK1/(2*M_PI);
71 for (i = 0; (i < n-nskip); i++)
75 lambda = eigval[i]*AMU;
76 w = sqrt(BOLTZMANN*temp/lambda)/NANO;
77 hwkT = (hbar*w)/(BOLTZMANN*temp);
78 dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
82 fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
83 i, w, lambda, hwkT, dS);
88 fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
92 fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
96 static void calc_entropy_schlitter(FILE *fp, int n, int nskip,
97 real *eigval, real temp)
103 double hbar, kt, kteh, S;
105 hbar = PLANCK1/(2*M_PI);
107 kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
110 fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
114 for (i = 0; (i < n-nskip); i++)
116 dd = 1+kteh*eigval[i];
121 fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S);
124 const char *proj_unit;
126 static real tick_spacing(real range, int minticks)
135 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
136 while (range/sp < minticks-1)
144 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
145 const char *title, const char *subtitle,
146 const char *xlabel, const char **ylabel,
147 int n, real *x, real **y, real ***sy,
148 real scale_x, gmx_bool bZero, gmx_bool bSplit,
149 const output_env_t oenv)
153 real min, max, xsp, ysp;
155 out = ffopen(file, "w");
156 if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
158 fprintf(out, "@ autoscale onread none\n");
160 for (g = 0; g < ngraphs; g++)
166 for (i = 0; i < n; i++)
182 for (s = 0; s < nsetspergraph; s++)
184 for (i = 0; i < n; i++)
186 if (sy[g][s][i] < min)
190 if (sy[g][s][i] > max)
203 min = min-0.1*(max-min);
205 max = max+0.1*(max-min);
206 xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
207 ysp = tick_spacing(max-min, 3);
208 if (output_env_get_print_xvgr_codes(oenv))
210 fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
213 fprintf(out, "@ title \"%s\"\n", title);
216 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
221 fprintf(out, "@ xaxis label \"%s\"\n", xlabel);
225 fprintf(out, "@ xaxis ticklabel off\n");
229 fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
230 fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
231 fprintf(out, "@ world ymin %g\n", min);
232 fprintf(out, "@ world ymax %g\n", max);
234 fprintf(out, "@ view xmin 0.15\n");
235 fprintf(out, "@ view xmax 0.85\n");
236 fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs);
237 fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs);
238 fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
239 fprintf(out, "@ xaxis tick major %g\n", xsp);
240 fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
241 fprintf(out, "@ xaxis ticklabel start type spec\n");
242 fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp);
243 fprintf(out, "@ yaxis tick major %g\n", ysp);
244 fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
245 fprintf(out, "@ yaxis ticklabel start type spec\n");
246 fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp);
247 if ((min < 0) && (max > 0))
249 fprintf(out, "@ zeroxaxis bar on\n");
250 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
253 for (s = 0; s < nsetspergraph; s++)
255 for (i = 0; i < n; i++)
257 if (bSplit && i > 0 && abs(x[i]) < 1e-5)
259 if (output_env_get_print_xvgr_codes(oenv))
264 fprintf(out, "%10.4f %10.5f\n",
265 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
267 if (output_env_get_print_xvgr_codes(oenv))
277 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
278 real *eigval1, int neig1, real *eigval2, int neig2)
282 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
286 n = min(n, min(neig1, neig2));
287 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
290 for (i = 0; i < n; i++)
297 eigval1[i] = sqrt(eigval1[i]);
300 for (i = n; i < neig1; i++)
302 trace1 += eigval1[i];
305 for (i = 0; i < n; i++)
312 eigval2[i] = sqrt(eigval2[i]);
315 for (i = n; i < neig2; i++)
317 trace2 += eigval2[i];
320 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
321 if (neig1 != n || neig2 != n)
323 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
324 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
326 fprintf(stdout, "Square root of the traces: %g and %g\n",
327 sqrt(sum1), sqrt(sum2));
330 for (i = 0; i < n; i++)
333 for (j = 0; j < n; j++)
336 for (k = 0; k < natoms; k++)
338 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
340 tmp += eigval2[j]*ip*ip;
342 sab += eigval1[i]*tmp;
345 samsb2 = sum1+sum2-2*sab;
351 fprintf(stdout, "The overlap of the covariance matrices:\n");
352 fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
353 tmp = 1-sab/sqrt(sum1*sum2);
358 fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
362 static void inprod_matrix(const char *matfile, int natoms,
363 int nvec1, int *eignr1, rvec **eigvec1,
364 int nvec2, int *eignr2, rvec **eigvec2,
365 gmx_bool bSelect, int noutvec, int *outvec)
369 int i, x1, y1, x, y, nlevels;
371 real inp, *t_x, *t_y, max;
379 for (y1 = 0; y1 < nx; y1++)
381 if (outvec[y1] < nvec2)
383 t_y[ny] = eignr2[outvec[y1]]+1;
392 for (y = 0; y < ny; y++)
394 t_y[y] = eignr2[y]+1;
398 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
404 for (x1 = 0; x1 < nx; x1++)
415 t_x[x1] = eignr1[x]+1;
416 fprintf(stderr, " %d", eignr1[x]+1);
417 for (y1 = 0; y1 < ny; y1++)
422 while (outvec[y1] >= nvec2)
432 for (i = 0; i < natoms; i++)
434 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
436 mat[x1][y1] = fabs(inp);
437 if (mat[x1][y1] > max)
443 fprintf(stderr, "\n");
444 rlo.r = 1; rlo.g = 1; rlo.b = 1;
445 rhi.r = 0; rhi.g = 0; rhi.b = 0;
447 out = ffopen(matfile, "w");
448 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
449 nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
453 static void overlap(const char *outfile, int natoms,
455 int nvec2, int *eignr2, rvec **eigvec2,
456 int noutvec, int *outvec,
457 const output_env_t oenv)
463 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
464 for (i = 0; i < noutvec; i++)
466 fprintf(stderr, "%d ", outvec[i]+1);
468 fprintf(stderr, "\n");
470 out = xvgropen(outfile, "Subspace overlap",
471 "Eigenvectors of trajectory 2", "Overlap", oenv);
472 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
475 for (x = 0; x < nvec2; x++)
477 for (v = 0; v < noutvec; v++)
481 for (i = 0; i < natoms; i++)
483 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
487 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
493 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
494 const char *projfile, const char *twodplotfile,
495 const char *threedplotfile, const char *filterfile, int skip,
496 const char *extremefile, gmx_bool bExtrAll, real extreme,
497 int nextr, t_atoms *atoms, int natoms, atom_id *index,
498 gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
499 real *sqrtm, rvec *xav,
500 int *eignr, rvec **eigvec,
501 int noutvec, int *outvec, gmx_bool bSplit,
502 const output_env_t oenv)
504 FILE *xvgrout = NULL;
505 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
506 t_trxstatus *out = NULL;
508 int noutvec_extr, imin, imax;
513 real t, inp, **inprod = NULL, min = 0, max = 0;
514 char str[STRLEN], str2[STRLEN], **ylabel, *c;
516 gmx_rmpbc_t gpbc = NULL;
522 noutvec_extr = noutvec;
532 snew(inprod, noutvec+1);
536 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
538 for (i = 0; i < noutvec; i++)
540 fprintf(stderr, "%d ", outvec[i]+1);
542 fprintf(stderr, "\n");
543 out = open_trx(filterfile, "w");
548 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
551 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);
557 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
560 for (i = 0; i < nat; i++)
570 gmx_rmpbc(gpbc, nat, box, xread);
572 if (nframes >= snew_size)
575 for (i = 0; i < noutvec+1; i++)
577 srenew(inprod[i], snew_size);
580 inprod[noutvec][nframes] = t;
581 /* calculate x: a fitted struture of the selected atoms */
584 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
585 do_fit(nat, w_rls, xref, xread);
587 for (i = 0; i < natoms; i++)
589 copy_rvec(xread[index[i]], x[i]);
592 for (v = 0; v < noutvec; v++)
595 /* calculate (mass-weighted) projection */
597 for (i = 0; i < natoms; i++)
599 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
600 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
601 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
603 inprod[v][nframes] = inp;
607 for (i = 0; i < natoms; i++)
609 for (d = 0; d < DIM; d++)
611 /* misuse xread for output */
612 xread[index[i]][d] = xav[i][d];
613 for (v = 0; v < noutvec; v++)
615 xread[index[i]][d] +=
616 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
620 write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
626 while (read_next_x(oenv, status, &t, xread, box));
636 snew(xread, atoms->nr);
641 gmx_rmpbc_done(gpbc);
647 snew(ylabel, noutvec);
648 for (v = 0; v < noutvec; v++)
650 sprintf(str, "vec %d", eignr[outvec[v]]+1);
651 ylabel[v] = strdup(str);
653 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
654 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
655 (const char **)ylabel,
656 nframes, inprod[noutvec], inprod, NULL,
657 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
662 sprintf(str, "projection on eigenvector %d (%s)",
663 eignr[outvec[0]]+1, proj_unit);
664 sprintf(str2, "projection on eigenvector %d (%s)",
665 eignr[outvec[noutvec-1]]+1, proj_unit);
666 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
668 for (i = 0; i < nframes; i++)
670 if (bSplit && i > 0 && abs(inprod[noutvec][i]) < 1e-5)
672 fprintf(xvgrout, "&\n");
674 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
685 char *resnm, *atnm, pdbform[STRLEN];
691 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
695 bPDB = fn2ftp(threedplotfile) == efPDB;
697 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
699 b4D = bPDB && (noutvec >= 4);
702 fprintf(stderr, "You have selected four or more eigenvectors:\n"
703 "fourth eigenvector will be plotted "
704 "in bfactor field of pdb file\n");
705 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
706 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
707 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
711 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
712 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
714 init_t_atoms(&atoms, nframes, FALSE);
718 resnm = strdup("PRJ");
722 fact = 10000.0/nframes;
729 for (i = 0; i < nframes; i++)
731 atoms.atomname[i] = &atnm;
732 atoms.atom[i].resind = i;
733 atoms.resinfo[i].name = &resnm;
734 atoms.resinfo[i].nr = ceil(i*fact);
735 atoms.resinfo[i].ic = ' ';
736 x[i][XX] = inprod[0][i];
737 x[i][YY] = inprod[1][i];
738 x[i][ZZ] = inprod[2][i];
744 if ( ( b4D || bSplit ) && bPDB)
746 strcpy(pdbform, get_pdbformat());
747 strcat(pdbform, "%8.4f%8.4f\n");
749 out = ffopen(threedplotfile, "w");
750 fprintf(out, "HEADER %s\n", str);
753 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
756 for (i = 0; i < atoms.nr; i++)
758 if (j > 0 && bSplit && abs(inprod[noutvec][i]) < 1e-5)
760 fprintf(out, "TER\n");
763 fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1,
764 PR_VEC(10*x[i]), 1.0, 10*b[i]);
767 fprintf(out, "CONECT%5d%5d\n", i, i+1);
771 fprintf(out, "TER\n");
776 write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
778 free_t_atoms(&atoms, FALSE);
783 snew(pmin, noutvec_extr);
784 snew(pmax, noutvec_extr);
787 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
789 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
792 for (v = 0; v < noutvec_extr; v++)
794 for (i = 0; i < nframes; i++)
796 if (inprod[v][i] < inprod[v][imin])
800 if (inprod[v][i] > inprod[v][imax])
805 pmin[v] = inprod[v][imin];
806 pmax[v] = inprod[v][imax];
807 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
809 pmin[v], imin, pmax[v], imax);
817 /* build format string for filename: */
818 strcpy(str, extremefile); /* copy filename */
819 c = strrchr(str, '.'); /* find where extention begins */
820 strcpy(str2, c); /* get extention */
821 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
822 for (v = 0; v < noutvec_extr; v++)
824 /* make filename using format string */
825 if (noutvec_extr == 1)
827 strcpy(str2, extremefile);
831 sprintf(str2, str, eignr[outvec[v]]+1);
833 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
834 nextr, outvec[v]+1, str2);
835 out = open_trx(str2, "w");
836 for (frame = 0; frame < nextr; frame++)
838 if ((extreme == 0) && (nextr <= 3))
840 for (i = 0; i < natoms; i++)
842 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
845 for (i = 0; i < natoms; i++)
847 for (d = 0; d < DIM; d++)
850 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
851 *eigvec[outvec[v]][i][d]/sqrtm[i]);
854 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
861 fprintf(stderr, "\n");
864 static void components(const char *outfile, int natoms,
865 int *eignr, rvec **eigvec,
866 int noutvec, int *outvec,
867 const output_env_t oenv)
871 char str[STRLEN], **ylabel;
873 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
875 snew(ylabel, noutvec);
878 for (i = 0; i < natoms; i++)
882 for (g = 0; g < noutvec; g++)
885 sprintf(str, "vec %d", eignr[v]+1);
886 ylabel[g] = strdup(str);
888 for (s = 0; s < 4; s++)
890 snew(y[g][s], natoms);
892 for (i = 0; i < natoms; i++)
894 y[g][0][i] = norm(eigvec[v][i]);
895 for (s = 0; s < 3; s++)
897 y[g][s+1][i] = eigvec[v][i][s];
901 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
902 "black: total, red: x, green: y, blue: z",
903 "Atom number", (const char **)ylabel,
904 natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
905 fprintf(stderr, "\n");
908 static void rmsf(const char *outfile, int natoms, real *sqrtm,
909 int *eignr, rvec **eigvec,
910 int noutvec, int *outvec,
911 real *eigval, int neig,
912 const output_env_t oenv)
916 char str[STRLEN], **ylabel;
918 for (i = 0; i < neig; i++)
926 fprintf(stderr, "Writing rmsf to %s\n", outfile);
928 snew(ylabel, noutvec);
931 for (i = 0; i < natoms; i++)
935 for (g = 0; g < noutvec; g++)
938 if (eignr[v] >= neig)
940 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
942 sprintf(str, "vec %d", eignr[v]+1);
943 ylabel[g] = strdup(str);
945 for (i = 0; i < natoms; i++)
947 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
950 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
951 "Atom number", (const char **)ylabel,
952 natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
953 fprintf(stderr, "\n");
956 int gmx_anaeig(int argc, char *argv[])
958 static const char *desc[] = {
959 "[TT]g_anaeig[tt] analyzes eigenvectors. The eigenvectors can be of a",
960 "covariance matrix ([TT]g_covar[tt]) or of a Normal Modes analysis",
961 "([TT]g_nmeig[tt]).[PAR]",
963 "When a trajectory is projected on eigenvectors, all structures are",
964 "fitted to the structure in the eigenvector file, if present, otherwise",
965 "to the structure in the structure file. When no run input file is",
966 "supplied, periodicity will not be taken into account. Most analyses",
967 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
968 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
970 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
971 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
973 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
974 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
976 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
977 "[TT]-first[tt] to [TT]-last[tt].",
978 "The projections of a trajectory on the eigenvectors of its",
979 "covariance matrix are called principal components (pc's).",
980 "It is often useful to check the cosine content of the pc's,",
981 "since the pc's of random diffusion are cosines with the number",
982 "of periods equal to half the pc index.",
983 "The cosine content of the pc's can be calculated with the program",
984 "[TT]g_analyze[tt].[PAR]",
986 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
987 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
989 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
990 "three selected eigenvectors.[PAR]",
992 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
993 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
995 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
996 "on the average structure and interpolate [TT]-nframes[tt] frames",
997 "between them, or set your own extremes with [TT]-max[tt]. The",
998 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
999 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1000 "will be written to separate files. Chain identifiers will be added",
1001 "when writing a [TT].pdb[tt] file with two or three structures (you",
1002 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1004 " Overlap calculations between covariance analysis:[BR]",
1005 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1007 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1008 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1009 "in file [TT]-v[tt].[PAR]",
1011 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1012 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1013 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1014 "have been set explicitly.[PAR]",
1016 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1017 "a single number for the overlap between the covariance matrices is",
1018 "generated. The formulas are:[BR]",
1019 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1020 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1021 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1022 "where M1 and M2 are the two covariance matrices and tr is the trace",
1023 "of a matrix. The numbers are proportional to the overlap of the square",
1024 "root of the fluctuations. The normalized overlap is the most useful",
1025 "number, it is 1 for identical matrices and 0 when the sampled",
1026 "subspaces are orthogonal.[PAR]",
1027 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1028 "computed based on the Quasiharmonic approach and based on",
1029 "Schlitter's formula."
1031 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1032 static real max = 0.0, temp = 298.15;
1033 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1035 { "-first", FALSE, etINT, {&first},
1036 "First eigenvector for analysis (-1 is select)" },
1037 { "-last", FALSE, etINT, {&last},
1038 "Last eigenvector for analysis (-1 is till the last)" },
1039 { "-skip", FALSE, etINT, {&skip},
1040 "Only analyse every nr-th frame" },
1041 { "-max", FALSE, etREAL, {&max},
1042 "Maximum for projection of the eigenvector on the average structure, "
1043 "max=0 gives the extremes" },
1044 { "-nframes", FALSE, etINT, {&nextr},
1045 "Number of frames for the extremes output" },
1046 { "-split", FALSE, etBOOL, {&bSplit},
1047 "Split eigenvector projections where time is zero" },
1048 { "-entropy", FALSE, etBOOL, {&bEntropy},
1049 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1050 { "-temp", FALSE, etREAL, {&temp},
1051 "Temperature for entropy calculations" },
1052 { "-nevskip", FALSE, etINT, {&nskip},
1053 "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." }
1055 #define NPA asize(pa)
1061 t_atoms *atoms = NULL;
1062 rvec *xtop, *xref1, *xref2, *xrefp = NULL;
1063 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1064 int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1065 rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1067 real xid, totmass, *sqrtm, *w_rls, t, lambda;
1070 const char *indexfile;
1073 int nout, *iout, noutvec, *outvec, nfit;
1074 atom_id *index, *ifit;
1075 const char *VecFile, *Vec2File, *topfile;
1076 const char *EigFile, *Eig2File;
1077 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1078 const char *TwoDPlotFile, *ThreeDPlotFile;
1079 const char *FilterFile, *ExtremeFile;
1080 const char *OverlapFile, *InpMatFile;
1081 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1082 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1083 real *eigval1 = NULL, *eigval2 = NULL;
1090 { efTRN, "-v", "eigenvec", ffREAD },
1091 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1092 { efTRX, "-f", NULL, ffOPTRD },
1093 { efTPS, NULL, NULL, ffOPTRD },
1094 { efNDX, NULL, NULL, ffOPTRD },
1095 { efXVG, "-eig", "eigenval", ffOPTRD },
1096 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1097 { efXVG, "-comp", "eigcomp", ffOPTWR },
1098 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1099 { efXVG, "-proj", "proj", ffOPTWR },
1100 { efXVG, "-2d", "2dproj", ffOPTWR },
1101 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1102 { efTRX, "-filt", "filtered", ffOPTWR },
1103 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1104 { efXVG, "-over", "overlap", ffOPTWR },
1105 { efXPM, "-inpr", "inprod", ffOPTWR }
1107 #define NFILE asize(fnm)
1109 parse_common_args(&argc, argv,
1110 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
1111 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv);
1113 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1115 VecFile = opt2fn("-v", NFILE, fnm);
1116 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1117 topfile = ftp2fn(efTPS, NFILE, fnm);
1118 EigFile = opt2fn_null("-eig", NFILE, fnm);
1119 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1120 CompFile = opt2fn_null("-comp", NFILE, fnm);
1121 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1122 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1123 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1124 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1125 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1126 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1127 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1128 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1130 bTop = fn2bTPX(topfile);
1131 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1132 || FilterFile || ExtremeFile;
1134 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1135 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1136 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1137 bVec2 = Vec2File || OverlapFile || InpMatFile;
1138 bM = RmsfFile || bProj;
1139 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1140 || TwoDPlotFile || ThreeDPlotFile;
1141 bIndex = bM || bProj;
1142 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1143 FilterFile || (bIndex && indexfile);
1144 bCompare = Vec2File || Eig2File;
1145 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1147 read_eigenvectors(VecFile, &natoms, &bFit1,
1148 &xref1, &bDMR1, &xav1, &bDMA1,
1149 &nvec1, &eignr1, &eigvec1, &eigval1);
1152 /* Overwrite eigenvalues from separate files if the user provides them */
1153 if (EigFile != NULL)
1155 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1156 if (neig_tmp != neig1)
1158 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1161 srenew(eigval1, neig1);
1162 for (j = 0; j < neig1; j++)
1164 real tmp = eigval1[j];
1165 eigval1[j] = xvgdata[1][j];
1166 if (debug && (eigval1[j] != tmp))
1168 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1169 j, tmp, eigval1[j]);
1172 for (j = 0; j < i; j++)
1177 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1184 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1186 calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1187 calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1194 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1196 read_eigenvectors(Vec2File, &neig2, &bFit2,
1197 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1202 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1206 if (Eig2File != NULL)
1208 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1209 srenew(eigval2, neig2);
1210 for (j = 0; j < neig2; j++)
1212 eigval2[j] = xvgdata[1][j];
1214 for (j = 0; j < i; j++)
1219 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1223 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1227 if ((xref1 == NULL) && (bM || bTraj))
1243 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1244 title, &top, &ePBC, &xtop, NULL, topbox, bM);
1246 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1247 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1248 /* Fitting is only required for the projection */
1253 printf("\nNote: the structure in %s should be the same\n"
1254 " as the one used for the fit in g_covar\n", topfile);
1256 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1257 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1259 snew(w_rls, atoms->nr);
1260 for (i = 0; (i < nfit); i++)
1264 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1268 w_rls[ifit[i]] = 1.0;
1272 snew(xrefp, atoms->nr);
1275 for (i = 0; (i < nfit); i++)
1277 copy_rvec(xref1[i], xrefp[ifit[i]]);
1282 /* The top coordinates are the fitting reference */
1283 for (i = 0; (i < nfit); i++)
1285 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1287 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1290 gmx_rmpbc_done(gpbc);
1295 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1296 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1299 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1304 snew(sqrtm, natoms);
1307 proj_unit = "u\\S1/2\\Nnm";
1308 for (i = 0; (i < natoms); i++)
1310 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1316 for (i = 0; (i < natoms); i++)
1326 for (i = 0; (i < natoms); i++)
1328 for (d = 0; (d < DIM); d++)
1330 t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1331 totmass += sqr(sqrtm[i]);
1334 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1335 " %.3f (nm)\n\n", sqrt(t/totmass));
1346 /* make an index from first to last */
1347 nout = last-first+1;
1349 for (i = 0; i < nout; i++)
1351 iout[i] = first-1+i;
1354 else if (ThreeDPlotFile)
1356 /* make an index of first+(0,1,2) and last */
1357 nout = bPDB3D ? 4 : 3;
1358 nout = min(last-first+1, nout);
1366 iout[nout-1] = last-1;
1370 /* make an index of first and last */
1379 printf("Select eigenvectors for output, end your selection with 0\n");
1386 srenew(iout, nout+1);
1387 if (1 != scanf("%d", &iout[nout]))
1389 gmx_fatal(FARGS, "Error reading user input");
1393 while (iout[nout] >= 0);
1397 /* make an index of the eigenvectors which are present */
1400 for (i = 0; i < nout; i++)
1403 while ((j < nvec1) && (eignr1[j] != iout[i]))
1407 if ((j < nvec1) && (eignr1[j] == iout[i]))
1409 outvec[noutvec] = j;
1413 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1416 fprintf(stderr, ":");
1417 for (j = 0; j < noutvec; j++)
1419 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1422 fprintf(stderr, "\n");
1426 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1431 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1437 project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1438 bTop ? &top : NULL, ePBC, topbox,
1439 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1440 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1441 bFit1, xrefp, nfit, ifit, w_rls,
1442 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1448 overlap(OverlapFile, natoms,
1449 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1454 inprod_matrix(InpMatFile, natoms,
1455 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1456 bFirstLastSet, noutvec, outvec);
1461 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1465 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1466 !bCompare && !bEntropy)
1468 fprintf(stderr, "\nIf you want some output,"
1469 " set one (or two or ...) of the output file options\n");
1473 view_all(oenv, NFILE, fnm);