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"
48 #include "gromacs/fileio/futil.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"
65 static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
68 double hwkT, w, dS, S = 0;
71 hbar = PLANCK1/(2*M_PI);
72 for (i = 0; (i < n-nskip); i++)
76 lambda = eigval[i]*AMU;
77 w = sqrt(BOLTZMANN*temp/lambda)/NANO;
78 hwkT = (hbar*w)/(BOLTZMANN*temp);
79 dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
83 fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
84 i, w, lambda, hwkT, dS);
89 fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
93 fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
97 static void calc_entropy_schlitter(FILE *fp, int n, int nskip,
98 real *eigval, real temp)
104 double hbar, kt, kteh, S;
106 hbar = PLANCK1/(2*M_PI);
108 kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
111 fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
115 for (i = 0; (i < n-nskip); i++)
117 dd = 1+kteh*eigval[i];
122 fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S);
125 const char *proj_unit;
127 static real tick_spacing(real range, int minticks)
136 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
137 while (range/sp < minticks-1)
145 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
146 const char *title, const char *subtitle,
147 const char *xlabel, const char **ylabel,
148 int n, real *x, real **y, real ***sy,
149 real scale_x, gmx_bool bZero, gmx_bool bSplit,
150 const output_env_t oenv)
154 real min, max, xsp, ysp;
156 out = ffopen(file, "w");
157 if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
159 fprintf(out, "@ autoscale onread none\n");
161 for (g = 0; g < ngraphs; g++)
167 for (i = 0; i < n; i++)
183 for (s = 0; s < nsetspergraph; s++)
185 for (i = 0; i < n; i++)
187 if (sy[g][s][i] < min)
191 if (sy[g][s][i] > max)
204 min = min-0.1*(max-min);
206 max = max+0.1*(max-min);
207 xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
208 ysp = tick_spacing(max-min, 3);
209 if (output_env_get_print_xvgr_codes(oenv))
211 fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
214 fprintf(out, "@ title \"%s\"\n", title);
217 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
222 fprintf(out, "@ xaxis label \"%s\"\n", xlabel);
226 fprintf(out, "@ xaxis ticklabel off\n");
230 fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
231 fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
232 fprintf(out, "@ world ymin %g\n", min);
233 fprintf(out, "@ world ymax %g\n", max);
235 fprintf(out, "@ view xmin 0.15\n");
236 fprintf(out, "@ view xmax 0.85\n");
237 fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs);
238 fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs);
239 fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
240 fprintf(out, "@ xaxis tick major %g\n", xsp);
241 fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
242 fprintf(out, "@ xaxis ticklabel start type spec\n");
243 fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp);
244 fprintf(out, "@ yaxis tick major %g\n", ysp);
245 fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
246 fprintf(out, "@ yaxis ticklabel start type spec\n");
247 fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp);
248 if ((min < 0) && (max > 0))
250 fprintf(out, "@ zeroxaxis bar on\n");
251 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
254 for (s = 0; s < nsetspergraph; s++)
256 for (i = 0; i < n; i++)
258 if (bSplit && i > 0 && abs(x[i]) < 1e-5)
260 if (output_env_get_print_xvgr_codes(oenv))
265 fprintf(out, "%10.4f %10.5f\n",
266 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
268 if (output_env_get_print_xvgr_codes(oenv))
278 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
279 real *eigval1, int neig1, real *eigval2, int neig2)
283 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
287 n = min(n, min(neig1, neig2));
288 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
291 for (i = 0; i < n; i++)
298 eigval1[i] = sqrt(eigval1[i]);
301 for (i = n; i < neig1; i++)
303 trace1 += eigval1[i];
306 for (i = 0; i < n; i++)
313 eigval2[i] = sqrt(eigval2[i]);
316 for (i = n; i < neig2; i++)
318 trace2 += eigval2[i];
321 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
322 if (neig1 != n || neig2 != n)
324 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
325 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
327 fprintf(stdout, "Square root of the traces: %g and %g\n",
328 sqrt(sum1), sqrt(sum2));
331 for (i = 0; i < n; i++)
334 for (j = 0; j < n; j++)
337 for (k = 0; k < natoms; k++)
339 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
341 tmp += eigval2[j]*ip*ip;
343 sab += eigval1[i]*tmp;
346 samsb2 = sum1+sum2-2*sab;
352 fprintf(stdout, "The overlap of the covariance matrices:\n");
353 fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
354 tmp = 1-sab/sqrt(sum1*sum2);
359 fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
363 static void inprod_matrix(const char *matfile, int natoms,
364 int nvec1, int *eignr1, rvec **eigvec1,
365 int nvec2, int *eignr2, rvec **eigvec2,
366 gmx_bool bSelect, int noutvec, int *outvec)
370 int i, x1, y1, x, y, nlevels;
372 real inp, *t_x, *t_y, max;
380 for (y1 = 0; y1 < nx; y1++)
382 if (outvec[y1] < nvec2)
384 t_y[ny] = eignr2[outvec[y1]]+1;
393 for (y = 0; y < ny; y++)
395 t_y[y] = eignr2[y]+1;
399 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
405 for (x1 = 0; x1 < nx; x1++)
416 t_x[x1] = eignr1[x]+1;
417 fprintf(stderr, " %d", eignr1[x]+1);
418 for (y1 = 0; y1 < ny; y1++)
423 while (outvec[y1] >= nvec2)
433 for (i = 0; i < natoms; i++)
435 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
437 mat[x1][y1] = fabs(inp);
438 if (mat[x1][y1] > max)
444 fprintf(stderr, "\n");
445 rlo.r = 1; rlo.g = 1; rlo.b = 1;
446 rhi.r = 0; rhi.g = 0; rhi.b = 0;
448 out = ffopen(matfile, "w");
449 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
450 nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
454 static void overlap(const char *outfile, int natoms,
456 int nvec2, int *eignr2, rvec **eigvec2,
457 int noutvec, int *outvec,
458 const output_env_t oenv)
464 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
465 for (i = 0; i < noutvec; i++)
467 fprintf(stderr, "%d ", outvec[i]+1);
469 fprintf(stderr, "\n");
471 out = xvgropen(outfile, "Subspace overlap",
472 "Eigenvectors of trajectory 2", "Overlap", oenv);
473 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
476 for (x = 0; x < nvec2; x++)
478 for (v = 0; v < noutvec; v++)
482 for (i = 0; i < natoms; i++)
484 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
488 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
494 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
495 const char *projfile, const char *twodplotfile,
496 const char *threedplotfile, const char *filterfile, int skip,
497 const char *extremefile, gmx_bool bExtrAll, real extreme,
498 int nextr, t_atoms *atoms, int natoms, atom_id *index,
499 gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
500 real *sqrtm, rvec *xav,
501 int *eignr, rvec **eigvec,
502 int noutvec, int *outvec, gmx_bool bSplit,
503 const output_env_t oenv)
505 FILE *xvgrout = NULL;
506 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
507 t_trxstatus *out = NULL;
509 int noutvec_extr, imin, imax;
514 real t, inp, **inprod = NULL, min = 0, max = 0;
515 char str[STRLEN], str2[STRLEN], **ylabel, *c;
517 gmx_rmpbc_t gpbc = NULL;
523 noutvec_extr = noutvec;
533 snew(inprod, noutvec+1);
537 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
539 for (i = 0; i < noutvec; i++)
541 fprintf(stderr, "%d ", outvec[i]+1);
543 fprintf(stderr, "\n");
544 out = open_trx(filterfile, "w");
549 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
552 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);
558 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
561 for (i = 0; i < nat; i++)
571 gmx_rmpbc(gpbc, nat, box, xread);
573 if (nframes >= snew_size)
576 for (i = 0; i < noutvec+1; i++)
578 srenew(inprod[i], snew_size);
581 inprod[noutvec][nframes] = t;
582 /* calculate x: a fitted struture of the selected atoms */
585 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
586 do_fit(nat, w_rls, xref, xread);
588 for (i = 0; i < natoms; i++)
590 copy_rvec(xread[index[i]], x[i]);
593 for (v = 0; v < noutvec; v++)
596 /* calculate (mass-weighted) projection */
598 for (i = 0; i < natoms; i++)
600 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
601 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
602 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
604 inprod[v][nframes] = inp;
608 for (i = 0; i < natoms; i++)
610 for (d = 0; d < DIM; d++)
612 /* misuse xread for output */
613 xread[index[i]][d] = xav[i][d];
614 for (v = 0; v < noutvec; v++)
616 xread[index[i]][d] +=
617 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
621 write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
627 while (read_next_x(oenv, status, &t, xread, box));
637 snew(xread, atoms->nr);
642 gmx_rmpbc_done(gpbc);
648 snew(ylabel, noutvec);
649 for (v = 0; v < noutvec; v++)
651 sprintf(str, "vec %d", eignr[outvec[v]]+1);
652 ylabel[v] = strdup(str);
654 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
655 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
656 (const char **)ylabel,
657 nframes, inprod[noutvec], inprod, NULL,
658 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
663 sprintf(str, "projection on eigenvector %d (%s)",
664 eignr[outvec[0]]+1, proj_unit);
665 sprintf(str2, "projection on eigenvector %d (%s)",
666 eignr[outvec[noutvec-1]]+1, proj_unit);
667 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
669 for (i = 0; i < nframes; i++)
671 if (bSplit && i > 0 && abs(inprod[noutvec][i]) < 1e-5)
673 fprintf(xvgrout, "&\n");
675 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
686 char *resnm, *atnm, pdbform[STRLEN];
692 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
696 bPDB = fn2ftp(threedplotfile) == efPDB;
698 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
700 b4D = bPDB && (noutvec >= 4);
703 fprintf(stderr, "You have selected four or more eigenvectors:\n"
704 "fourth eigenvector will be plotted "
705 "in bfactor field of pdb file\n");
706 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
707 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
708 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
712 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
713 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
715 init_t_atoms(&atoms, nframes, FALSE);
719 resnm = strdup("PRJ");
723 fact = 10000.0/nframes;
730 for (i = 0; i < nframes; i++)
732 atoms.atomname[i] = &atnm;
733 atoms.atom[i].resind = i;
734 atoms.resinfo[i].name = &resnm;
735 atoms.resinfo[i].nr = ceil(i*fact);
736 atoms.resinfo[i].ic = ' ';
737 x[i][XX] = inprod[0][i];
738 x[i][YY] = inprod[1][i];
739 x[i][ZZ] = inprod[2][i];
745 if ( ( b4D || bSplit ) && bPDB)
747 strcpy(pdbform, get_pdbformat());
748 strcat(pdbform, "%8.4f%8.4f\n");
750 out = ffopen(threedplotfile, "w");
751 fprintf(out, "HEADER %s\n", str);
754 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
757 for (i = 0; i < atoms.nr; i++)
759 if (j > 0 && bSplit && abs(inprod[noutvec][i]) < 1e-5)
761 fprintf(out, "TER\n");
764 fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1,
765 PR_VEC(10*x[i]), 1.0, 10*b[i]);
768 fprintf(out, "CONECT%5d%5d\n", i, i+1);
772 fprintf(out, "TER\n");
777 write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
779 free_t_atoms(&atoms, FALSE);
784 snew(pmin, noutvec_extr);
785 snew(pmax, noutvec_extr);
788 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
790 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
793 for (v = 0; v < noutvec_extr; v++)
795 for (i = 0; i < nframes; i++)
797 if (inprod[v][i] < inprod[v][imin])
801 if (inprod[v][i] > inprod[v][imax])
806 pmin[v] = inprod[v][imin];
807 pmax[v] = inprod[v][imax];
808 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
810 pmin[v], imin, pmax[v], imax);
818 /* build format string for filename: */
819 strcpy(str, extremefile); /* copy filename */
820 c = strrchr(str, '.'); /* find where extention begins */
821 strcpy(str2, c); /* get extention */
822 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
823 for (v = 0; v < noutvec_extr; v++)
825 /* make filename using format string */
826 if (noutvec_extr == 1)
828 strcpy(str2, extremefile);
832 sprintf(str2, str, eignr[outvec[v]]+1);
834 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
835 nextr, outvec[v]+1, str2);
836 out = open_trx(str2, "w");
837 for (frame = 0; frame < nextr; frame++)
839 if ((extreme == 0) && (nextr <= 3))
841 for (i = 0; i < natoms; i++)
843 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
846 for (i = 0; i < natoms; i++)
848 for (d = 0; d < DIM; d++)
851 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
852 *eigvec[outvec[v]][i][d]/sqrtm[i]);
855 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
862 fprintf(stderr, "\n");
865 static void components(const char *outfile, int natoms,
866 int *eignr, rvec **eigvec,
867 int noutvec, int *outvec,
868 const output_env_t oenv)
872 char str[STRLEN], **ylabel;
874 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
876 snew(ylabel, noutvec);
879 for (i = 0; i < natoms; i++)
883 for (g = 0; g < noutvec; g++)
886 sprintf(str, "vec %d", eignr[v]+1);
887 ylabel[g] = strdup(str);
889 for (s = 0; s < 4; s++)
891 snew(y[g][s], natoms);
893 for (i = 0; i < natoms; i++)
895 y[g][0][i] = norm(eigvec[v][i]);
896 for (s = 0; s < 3; s++)
898 y[g][s+1][i] = eigvec[v][i][s];
902 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
903 "black: total, red: x, green: y, blue: z",
904 "Atom number", (const char **)ylabel,
905 natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
906 fprintf(stderr, "\n");
909 static void rmsf(const char *outfile, int natoms, real *sqrtm,
910 int *eignr, rvec **eigvec,
911 int noutvec, int *outvec,
912 real *eigval, int neig,
913 const output_env_t oenv)
917 char str[STRLEN], **ylabel;
919 for (i = 0; i < neig; i++)
927 fprintf(stderr, "Writing rmsf to %s\n", outfile);
929 snew(ylabel, noutvec);
932 for (i = 0; i < natoms; i++)
936 for (g = 0; g < noutvec; g++)
939 if (eignr[v] >= neig)
941 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
943 sprintf(str, "vec %d", eignr[v]+1);
944 ylabel[g] = strdup(str);
946 for (i = 0; i < natoms; i++)
948 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
951 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
952 "Atom number", (const char **)ylabel,
953 natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
954 fprintf(stderr, "\n");
957 int gmx_anaeig(int argc, char *argv[])
959 static const char *desc[] = {
960 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
961 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
962 "([gmx-nmeig]).[PAR]",
964 "When a trajectory is projected on eigenvectors, all structures are",
965 "fitted to the structure in the eigenvector file, if present, otherwise",
966 "to the structure in the structure file. When no run input file is",
967 "supplied, periodicity will not be taken into account. Most analyses",
968 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
969 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
971 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
972 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
974 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
975 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
977 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
978 "[TT]-first[tt] to [TT]-last[tt].",
979 "The projections of a trajectory on the eigenvectors of its",
980 "covariance matrix are called principal components (pc's).",
981 "It is often useful to check the cosine content of the pc's,",
982 "since the pc's of random diffusion are cosines with the number",
983 "of periods equal to half the pc index.",
984 "The cosine content of the pc's can be calculated with the program",
985 "[gmx-analyze].[PAR]",
987 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
988 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
990 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
991 "three selected eigenvectors.[PAR]",
993 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
994 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
996 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
997 "on the average structure and interpolate [TT]-nframes[tt] frames",
998 "between them, or set your own extremes with [TT]-max[tt]. The",
999 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
1000 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1001 "will be written to separate files. Chain identifiers will be added",
1002 "when writing a [TT].pdb[tt] file with two or three structures (you",
1003 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1005 " Overlap calculations between covariance analysis:[BR]",
1006 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1008 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1009 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1010 "in file [TT]-v[tt].[PAR]",
1012 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1013 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1014 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1015 "have been set explicitly.[PAR]",
1017 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1018 "a single number for the overlap between the covariance matrices is",
1019 "generated. The formulas are:[BR]",
1020 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1021 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1022 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1023 "where M1 and M2 are the two covariance matrices and tr is the trace",
1024 "of a matrix. The numbers are proportional to the overlap of the square",
1025 "root of the fluctuations. The normalized overlap is the most useful",
1026 "number, it is 1 for identical matrices and 0 when the sampled",
1027 "subspaces are orthogonal.[PAR]",
1028 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1029 "computed based on the Quasiharmonic approach and based on",
1030 "Schlitter's formula."
1032 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1033 static real max = 0.0, temp = 298.15;
1034 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1036 { "-first", FALSE, etINT, {&first},
1037 "First eigenvector for analysis (-1 is select)" },
1038 { "-last", FALSE, etINT, {&last},
1039 "Last eigenvector for analysis (-1 is till the last)" },
1040 { "-skip", FALSE, etINT, {&skip},
1041 "Only analyse every nr-th frame" },
1042 { "-max", FALSE, etREAL, {&max},
1043 "Maximum for projection of the eigenvector on the average structure, "
1044 "max=0 gives the extremes" },
1045 { "-nframes", FALSE, etINT, {&nextr},
1046 "Number of frames for the extremes output" },
1047 { "-split", FALSE, etBOOL, {&bSplit},
1048 "Split eigenvector projections where time is zero" },
1049 { "-entropy", FALSE, etBOOL, {&bEntropy},
1050 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1051 { "-temp", FALSE, etREAL, {&temp},
1052 "Temperature for entropy calculations" },
1053 { "-nevskip", FALSE, etINT, {&nskip},
1054 "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." }
1056 #define NPA asize(pa)
1062 t_atoms *atoms = NULL;
1063 rvec *xtop, *xref1, *xref2, *xrefp = NULL;
1064 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1065 int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1066 rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1068 real xid, totmass, *sqrtm, *w_rls, t, lambda;
1071 const char *indexfile;
1074 int nout, *iout, noutvec, *outvec, nfit;
1075 atom_id *index, *ifit;
1076 const char *VecFile, *Vec2File, *topfile;
1077 const char *EigFile, *Eig2File;
1078 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1079 const char *TwoDPlotFile, *ThreeDPlotFile;
1080 const char *FilterFile, *ExtremeFile;
1081 const char *OverlapFile, *InpMatFile;
1082 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1083 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1084 real *eigval1 = NULL, *eigval2 = NULL;
1091 { efTRN, "-v", "eigenvec", ffREAD },
1092 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1093 { efTRX, "-f", NULL, ffOPTRD },
1094 { efTPS, NULL, NULL, ffOPTRD },
1095 { efNDX, NULL, NULL, ffOPTRD },
1096 { efXVG, "-eig", "eigenval", ffOPTRD },
1097 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1098 { efXVG, "-comp", "eigcomp", ffOPTWR },
1099 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1100 { efXVG, "-proj", "proj", ffOPTWR },
1101 { efXVG, "-2d", "2dproj", ffOPTWR },
1102 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1103 { efTRX, "-filt", "filtered", ffOPTWR },
1104 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1105 { efXVG, "-over", "overlap", ffOPTWR },
1106 { efXPM, "-inpr", "inprod", ffOPTWR }
1108 #define NFILE asize(fnm)
1110 if (!parse_common_args(&argc, argv,
1111 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
1112 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1117 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1119 VecFile = opt2fn("-v", NFILE, fnm);
1120 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1121 topfile = ftp2fn(efTPS, NFILE, fnm);
1122 EigFile = opt2fn_null("-eig", NFILE, fnm);
1123 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1124 CompFile = opt2fn_null("-comp", NFILE, fnm);
1125 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1126 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1127 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1128 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1129 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1130 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1131 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1132 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1134 bTop = fn2bTPX(topfile);
1135 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1136 || FilterFile || ExtremeFile;
1138 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1139 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1140 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1141 bVec2 = Vec2File || OverlapFile || InpMatFile;
1142 bM = RmsfFile || bProj;
1143 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1144 || TwoDPlotFile || ThreeDPlotFile;
1145 bIndex = bM || bProj;
1146 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1147 FilterFile || (bIndex && indexfile);
1148 bCompare = Vec2File || Eig2File;
1149 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1151 read_eigenvectors(VecFile, &natoms, &bFit1,
1152 &xref1, &bDMR1, &xav1, &bDMA1,
1153 &nvec1, &eignr1, &eigvec1, &eigval1);
1156 /* Overwrite eigenvalues from separate files if the user provides them */
1157 if (EigFile != NULL)
1159 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1160 if (neig_tmp != neig1)
1162 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1165 srenew(eigval1, neig1);
1166 for (j = 0; j < neig1; j++)
1168 real tmp = eigval1[j];
1169 eigval1[j] = xvgdata[1][j];
1170 if (debug && (eigval1[j] != tmp))
1172 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1173 j, tmp, eigval1[j]);
1176 for (j = 0; j < i; j++)
1181 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1188 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1190 calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1191 calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1198 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1200 read_eigenvectors(Vec2File, &neig2, &bFit2,
1201 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1206 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1210 if (Eig2File != NULL)
1212 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1213 srenew(eigval2, neig2);
1214 for (j = 0; j < neig2; j++)
1216 eigval2[j] = xvgdata[1][j];
1218 for (j = 0; j < i; j++)
1223 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1227 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1231 if ((xref1 == NULL) && (bM || bTraj))
1247 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1248 title, &top, &ePBC, &xtop, NULL, topbox, bM);
1250 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1251 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1252 /* Fitting is only required for the projection */
1257 printf("\nNote: the structure in %s should be the same\n"
1258 " as the one used for the fit in g_covar\n", topfile);
1260 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1261 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1263 snew(w_rls, atoms->nr);
1264 for (i = 0; (i < nfit); i++)
1268 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1272 w_rls[ifit[i]] = 1.0;
1276 snew(xrefp, atoms->nr);
1279 for (i = 0; (i < nfit); i++)
1281 copy_rvec(xref1[i], xrefp[ifit[i]]);
1286 /* The top coordinates are the fitting reference */
1287 for (i = 0; (i < nfit); i++)
1289 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1291 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1294 gmx_rmpbc_done(gpbc);
1299 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1300 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1303 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1308 snew(sqrtm, natoms);
1311 proj_unit = "u\\S1/2\\Nnm";
1312 for (i = 0; (i < natoms); i++)
1314 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1320 for (i = 0; (i < natoms); i++)
1330 for (i = 0; (i < natoms); i++)
1332 for (d = 0; (d < DIM); d++)
1334 t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1335 totmass += sqr(sqrtm[i]);
1338 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1339 " %.3f (nm)\n\n", sqrt(t/totmass));
1350 /* make an index from first to last */
1351 nout = last-first+1;
1353 for (i = 0; i < nout; i++)
1355 iout[i] = first-1+i;
1358 else if (ThreeDPlotFile)
1360 /* make an index of first+(0,1,2) and last */
1361 nout = bPDB3D ? 4 : 3;
1362 nout = min(last-first+1, nout);
1370 iout[nout-1] = last-1;
1374 /* make an index of first and last */
1383 printf("Select eigenvectors for output, end your selection with 0\n");
1390 srenew(iout, nout+1);
1391 if (1 != scanf("%d", &iout[nout]))
1393 gmx_fatal(FARGS, "Error reading user input");
1397 while (iout[nout] >= 0);
1401 /* make an index of the eigenvectors which are present */
1404 for (i = 0; i < nout; i++)
1407 while ((j < nvec1) && (eignr1[j] != iout[i]))
1411 if ((j < nvec1) && (eignr1[j] == iout[i]))
1413 outvec[noutvec] = j;
1417 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1420 fprintf(stderr, ":");
1421 for (j = 0; j < noutvec; j++)
1423 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1426 fprintf(stderr, "\n");
1430 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1435 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1441 project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1442 bTop ? &top : NULL, ePBC, topbox,
1443 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1444 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1445 bFit1, xrefp, nfit, ifit, w_rls,
1446 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1452 overlap(OverlapFile, natoms,
1453 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1458 inprod_matrix(InpMatFile, natoms,
1459 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1460 bFirstLastSet, noutvec, outvec);
1465 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1469 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1470 !bCompare && !bEntropy)
1472 fprintf(stderr, "\nIf you want some output,"
1473 " set one (or two or ...) of the output file options\n");
1477 view_all(oenv, NFILE, fnm);