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.
42 #include "gromacs/commandline/pargs.h"
47 #include "gmx_fatal.h"
50 #include "gromacs/fileio/futil.h"
52 #include "gromacs/fileio/pdbio.h"
53 #include "gromacs/fileio/confio.h"
54 #include "gromacs/fileio/tpxio.h"
55 #include "gromacs/fileio/trxio.h"
56 #include "gromacs/fileio/matio.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 && abs(x[i]) < 1e-5)
261 if (output_env_get_print_xvgr_codes(oenv))
266 fprintf(out, "%10.4f %10.5f\n",
267 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
269 if (output_env_get_print_xvgr_codes(oenv))
279 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
280 real *eigval1, int neig1, real *eigval2, int neig2)
284 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
288 n = min(n, min(neig1, neig2));
289 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
292 for (i = 0; i < n; i++)
299 eigval1[i] = sqrt(eigval1[i]);
302 for (i = n; i < neig1; i++)
304 trace1 += eigval1[i];
307 for (i = 0; i < n; i++)
314 eigval2[i] = sqrt(eigval2[i]);
317 for (i = n; i < neig2; i++)
319 trace2 += eigval2[i];
322 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
323 if (neig1 != n || neig2 != n)
325 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
326 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
328 fprintf(stdout, "Square root of the traces: %g and %g\n",
329 sqrt(sum1), sqrt(sum2));
332 for (i = 0; i < n; i++)
335 for (j = 0; j < n; j++)
338 for (k = 0; k < natoms; k++)
340 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
342 tmp += eigval2[j]*ip*ip;
344 sab += eigval1[i]*tmp;
347 samsb2 = sum1+sum2-2*sab;
353 fprintf(stdout, "The overlap of the covariance matrices:\n");
354 fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
355 tmp = 1-sab/sqrt(sum1*sum2);
360 fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
364 static void inprod_matrix(const char *matfile, int natoms,
365 int nvec1, int *eignr1, rvec **eigvec1,
366 int nvec2, int *eignr2, rvec **eigvec2,
367 gmx_bool bSelect, int noutvec, int *outvec)
371 int i, x1, y1, x, y, nlevels;
373 real inp, *t_x, *t_y, max;
381 for (y1 = 0; y1 < nx; y1++)
383 if (outvec[y1] < nvec2)
385 t_y[ny] = eignr2[outvec[y1]]+1;
394 for (y = 0; y < ny; y++)
396 t_y[y] = eignr2[y]+1;
400 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
406 for (x1 = 0; x1 < nx; x1++)
417 t_x[x1] = eignr1[x]+1;
418 fprintf(stderr, " %d", eignr1[x]+1);
419 for (y1 = 0; y1 < ny; y1++)
424 while (outvec[y1] >= nvec2)
434 for (i = 0; i < natoms; i++)
436 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
438 mat[x1][y1] = fabs(inp);
439 if (mat[x1][y1] > max)
445 fprintf(stderr, "\n");
446 rlo.r = 1; rlo.g = 1; rlo.b = 1;
447 rhi.r = 0; rhi.g = 0; rhi.b = 0;
449 out = gmx_ffopen(matfile, "w");
450 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
451 nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
455 static void overlap(const char *outfile, int natoms,
457 int nvec2, int *eignr2, rvec **eigvec2,
458 int noutvec, int *outvec,
459 const output_env_t oenv)
465 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
466 for (i = 0; i < noutvec; i++)
468 fprintf(stderr, "%d ", outvec[i]+1);
470 fprintf(stderr, "\n");
472 out = xvgropen(outfile, "Subspace overlap",
473 "Eigenvectors of trajectory 2", "Overlap", oenv);
474 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
477 for (x = 0; x < nvec2; x++)
479 for (v = 0; v < noutvec; v++)
483 for (i = 0; i < natoms; i++)
485 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
489 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
495 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
496 const char *projfile, const char *twodplotfile,
497 const char *threedplotfile, const char *filterfile, int skip,
498 const char *extremefile, gmx_bool bExtrAll, real extreme,
499 int nextr, t_atoms *atoms, int natoms, atom_id *index,
500 gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
501 real *sqrtm, rvec *xav,
502 int *eignr, rvec **eigvec,
503 int noutvec, int *outvec, gmx_bool bSplit,
504 const output_env_t oenv)
506 FILE *xvgrout = NULL;
507 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
508 t_trxstatus *out = NULL;
510 int noutvec_extr, imin, imax;
515 real t, inp, **inprod = NULL, min = 0, max = 0;
516 char str[STRLEN], str2[STRLEN], **ylabel, *c;
518 gmx_rmpbc_t gpbc = NULL;
524 noutvec_extr = noutvec;
534 snew(inprod, noutvec+1);
538 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
540 for (i = 0; i < noutvec; i++)
542 fprintf(stderr, "%d ", outvec[i]+1);
544 fprintf(stderr, "\n");
545 out = open_trx(filterfile, "w");
550 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
553 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);
559 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
562 for (i = 0; i < nat; i++)
572 gmx_rmpbc(gpbc, nat, box, xread);
574 if (nframes >= snew_size)
577 for (i = 0; i < noutvec+1; i++)
579 srenew(inprod[i], snew_size);
582 inprod[noutvec][nframes] = t;
583 /* calculate x: a fitted struture of the selected atoms */
586 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
587 do_fit(nat, w_rls, xref, xread);
589 for (i = 0; i < natoms; i++)
591 copy_rvec(xread[index[i]], x[i]);
594 for (v = 0; v < noutvec; v++)
597 /* calculate (mass-weighted) projection */
599 for (i = 0; i < natoms; i++)
601 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
602 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
603 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
605 inprod[v][nframes] = inp;
609 for (i = 0; i < natoms; i++)
611 for (d = 0; d < DIM; d++)
613 /* misuse xread for output */
614 xread[index[i]][d] = xav[i][d];
615 for (v = 0; v < noutvec; v++)
617 xread[index[i]][d] +=
618 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
622 write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
628 while (read_next_x(oenv, status, &t, xread, box));
638 snew(xread, atoms->nr);
643 gmx_rmpbc_done(gpbc);
649 snew(ylabel, noutvec);
650 for (v = 0; v < noutvec; v++)
652 sprintf(str, "vec %d", eignr[outvec[v]]+1);
653 ylabel[v] = strdup(str);
655 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
656 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
657 (const char **)ylabel,
658 nframes, inprod[noutvec], inprod, NULL,
659 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
664 sprintf(str, "projection on eigenvector %d (%s)",
665 eignr[outvec[0]]+1, proj_unit);
666 sprintf(str2, "projection on eigenvector %d (%s)",
667 eignr[outvec[noutvec-1]]+1, proj_unit);
668 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
670 for (i = 0; i < nframes; i++)
672 if (bSplit && i > 0 && abs(inprod[noutvec][i]) < 1e-5)
674 fprintf(xvgrout, "&\n");
676 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
678 gmx_ffclose(xvgrout);
687 char *resnm, *atnm, pdbform[STRLEN];
693 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
697 bPDB = fn2ftp(threedplotfile) == efPDB;
699 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
701 b4D = bPDB && (noutvec >= 4);
704 fprintf(stderr, "You have selected four or more eigenvectors:\n"
705 "fourth eigenvector will be plotted "
706 "in bfactor field of pdb file\n");
707 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
708 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
709 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
713 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
714 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
716 init_t_atoms(&atoms, nframes, FALSE);
720 resnm = strdup("PRJ");
724 fact = 10000.0/nframes;
731 for (i = 0; i < nframes; i++)
733 atoms.atomname[i] = &atnm;
734 atoms.atom[i].resind = i;
735 atoms.resinfo[i].name = &resnm;
736 atoms.resinfo[i].nr = ceil(i*fact);
737 atoms.resinfo[i].ic = ' ';
738 x[i][XX] = inprod[0][i];
739 x[i][YY] = inprod[1][i];
740 x[i][ZZ] = inprod[2][i];
746 if ( ( b4D || bSplit ) && bPDB)
748 strcpy(pdbform, get_pdbformat());
749 strcat(pdbform, "%8.4f%8.4f\n");
751 out = gmx_ffopen(threedplotfile, "w");
752 fprintf(out, "HEADER %s\n", str);
755 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
758 for (i = 0; i < atoms.nr; i++)
760 if (j > 0 && bSplit && abs(inprod[noutvec][i]) < 1e-5)
762 fprintf(out, "TER\n");
765 fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1,
766 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i]);
769 fprintf(out, "CONECT%5d%5d\n", i, i+1);
773 fprintf(out, "TER\n");
778 write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
780 free_t_atoms(&atoms, FALSE);
785 snew(pmin, noutvec_extr);
786 snew(pmax, noutvec_extr);
789 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
791 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
794 for (v = 0; v < noutvec_extr; v++)
796 for (i = 0; i < nframes; i++)
798 if (inprod[v][i] < inprod[v][imin])
802 if (inprod[v][i] > inprod[v][imax])
807 pmin[v] = inprod[v][imin];
808 pmax[v] = inprod[v][imax];
809 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
811 pmin[v], imin, pmax[v], imax);
819 /* build format string for filename: */
820 strcpy(str, extremefile); /* copy filename */
821 c = strrchr(str, '.'); /* find where extention begins */
822 strcpy(str2, c); /* get extention */
823 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
824 for (v = 0; v < noutvec_extr; v++)
826 /* make filename using format string */
827 if (noutvec_extr == 1)
829 strcpy(str2, extremefile);
833 sprintf(str2, str, eignr[outvec[v]]+1);
835 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
836 nextr, outvec[v]+1, str2);
837 out = open_trx(str2, "w");
838 for (frame = 0; frame < nextr; frame++)
840 if ((extreme == 0) && (nextr <= 3))
842 for (i = 0; i < natoms; i++)
844 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
847 for (i = 0; i < natoms; i++)
849 for (d = 0; d < DIM; d++)
852 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
853 *eigvec[outvec[v]][i][d]/sqrtm[i]);
856 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
863 fprintf(stderr, "\n");
866 static void components(const char *outfile, int natoms,
867 int *eignr, rvec **eigvec,
868 int noutvec, int *outvec,
869 const output_env_t oenv)
873 char str[STRLEN], **ylabel;
875 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
877 snew(ylabel, noutvec);
880 for (i = 0; i < natoms; i++)
884 for (g = 0; g < noutvec; g++)
887 sprintf(str, "vec %d", eignr[v]+1);
888 ylabel[g] = strdup(str);
890 for (s = 0; s < 4; s++)
892 snew(y[g][s], natoms);
894 for (i = 0; i < natoms; i++)
896 y[g][0][i] = norm(eigvec[v][i]);
897 for (s = 0; s < 3; s++)
899 y[g][s+1][i] = eigvec[v][i][s];
903 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
904 "black: total, red: x, green: y, blue: z",
905 "Atom number", (const char **)ylabel,
906 natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
907 fprintf(stderr, "\n");
910 static void rmsf(const char *outfile, int natoms, real *sqrtm,
911 int *eignr, rvec **eigvec,
912 int noutvec, int *outvec,
913 real *eigval, int neig,
914 const output_env_t oenv)
918 char str[STRLEN], **ylabel;
920 for (i = 0; i < neig; i++)
928 fprintf(stderr, "Writing rmsf to %s\n", outfile);
930 snew(ylabel, noutvec);
933 for (i = 0; i < natoms; i++)
937 for (g = 0; g < noutvec; g++)
940 if (eignr[v] >= neig)
942 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
944 sprintf(str, "vec %d", eignr[v]+1);
945 ylabel[g] = strdup(str);
947 for (i = 0; i < natoms; i++)
949 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
952 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
953 "Atom number", (const char **)ylabel,
954 natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
955 fprintf(stderr, "\n");
958 int gmx_anaeig(int argc, char *argv[])
960 static const char *desc[] = {
961 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
962 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
963 "([gmx-nmeig]).[PAR]",
965 "When a trajectory is projected on eigenvectors, all structures are",
966 "fitted to the structure in the eigenvector file, if present, otherwise",
967 "to the structure in the structure file. When no run input file is",
968 "supplied, periodicity will not be taken into account. Most analyses",
969 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
970 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
972 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
973 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
975 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
976 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
978 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
979 "[TT]-first[tt] to [TT]-last[tt].",
980 "The projections of a trajectory on the eigenvectors of its",
981 "covariance matrix are called principal components (pc's).",
982 "It is often useful to check the cosine content of the pc's,",
983 "since the pc's of random diffusion are cosines with the number",
984 "of periods equal to half the pc index.",
985 "The cosine content of the pc's can be calculated with the program",
986 "[gmx-analyze].[PAR]",
988 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
989 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
991 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
992 "three selected eigenvectors.[PAR]",
994 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
995 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
997 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
998 "on the average structure and interpolate [TT]-nframes[tt] frames",
999 "between them, or set your own extremes with [TT]-max[tt]. The",
1000 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
1001 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1002 "will be written to separate files. Chain identifiers will be added",
1003 "when writing a [TT].pdb[tt] file with two or three structures (you",
1004 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1006 " Overlap calculations between covariance analysis:[BR]",
1007 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1009 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1010 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1011 "in file [TT]-v[tt].[PAR]",
1013 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1014 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1015 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1016 "have been set explicitly.[PAR]",
1018 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1019 "a single number for the overlap between the covariance matrices is",
1020 "generated. The formulas are:[BR]",
1021 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1022 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1023 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1024 "where M1 and M2 are the two covariance matrices and tr is the trace",
1025 "of a matrix. The numbers are proportional to the overlap of the square",
1026 "root of the fluctuations. The normalized overlap is the most useful",
1027 "number, it is 1 for identical matrices and 0 when the sampled",
1028 "subspaces are orthogonal.[PAR]",
1029 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1030 "computed based on the Quasiharmonic approach and based on",
1031 "Schlitter's formula."
1033 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1034 static real max = 0.0, temp = 298.15;
1035 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1037 { "-first", FALSE, etINT, {&first},
1038 "First eigenvector for analysis (-1 is select)" },
1039 { "-last", FALSE, etINT, {&last},
1040 "Last eigenvector for analysis (-1 is till the last)" },
1041 { "-skip", FALSE, etINT, {&skip},
1042 "Only analyse every nr-th frame" },
1043 { "-max", FALSE, etREAL, {&max},
1044 "Maximum for projection of the eigenvector on the average structure, "
1045 "max=0 gives the extremes" },
1046 { "-nframes", FALSE, etINT, {&nextr},
1047 "Number of frames for the extremes output" },
1048 { "-split", FALSE, etBOOL, {&bSplit},
1049 "Split eigenvector projections where time is zero" },
1050 { "-entropy", FALSE, etBOOL, {&bEntropy},
1051 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1052 { "-temp", FALSE, etREAL, {&temp},
1053 "Temperature for entropy calculations" },
1054 { "-nevskip", FALSE, etINT, {&nskip},
1055 "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." }
1057 #define NPA asize(pa)
1063 t_atoms *atoms = NULL;
1064 rvec *xtop, *xref1, *xref2, *xrefp = NULL;
1065 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1066 int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1067 rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1069 real xid, totmass, *sqrtm, *w_rls, t, lambda;
1072 const char *indexfile;
1075 int nout, *iout, noutvec, *outvec, nfit;
1076 atom_id *index, *ifit;
1077 const char *VecFile, *Vec2File, *topfile;
1078 const char *EigFile, *Eig2File;
1079 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1080 const char *TwoDPlotFile, *ThreeDPlotFile;
1081 const char *FilterFile, *ExtremeFile;
1082 const char *OverlapFile, *InpMatFile;
1083 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1084 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1085 real *eigval1 = NULL, *eigval2 = NULL;
1092 { efTRN, "-v", "eigenvec", ffREAD },
1093 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1094 { efTRX, "-f", NULL, ffOPTRD },
1095 { efTPS, NULL, NULL, ffOPTRD },
1096 { efNDX, NULL, NULL, ffOPTRD },
1097 { efXVG, "-eig", "eigenval", ffOPTRD },
1098 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1099 { efXVG, "-comp", "eigcomp", ffOPTWR },
1100 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1101 { efXVG, "-proj", "proj", ffOPTWR },
1102 { efXVG, "-2d", "2dproj", ffOPTWR },
1103 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1104 { efTRX, "-filt", "filtered", ffOPTWR },
1105 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1106 { efXVG, "-over", "overlap", ffOPTWR },
1107 { efXPM, "-inpr", "inprod", ffOPTWR }
1109 #define NFILE asize(fnm)
1111 if (!parse_common_args(&argc, argv,
1112 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
1113 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1118 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1120 VecFile = opt2fn("-v", NFILE, fnm);
1121 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1122 topfile = ftp2fn(efTPS, NFILE, fnm);
1123 EigFile = opt2fn_null("-eig", NFILE, fnm);
1124 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1125 CompFile = opt2fn_null("-comp", NFILE, fnm);
1126 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1127 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1128 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1129 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1130 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1131 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1132 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1133 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1135 bTop = fn2bTPX(topfile);
1136 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1137 || FilterFile || ExtremeFile;
1139 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1140 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1141 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1142 bVec2 = Vec2File || OverlapFile || InpMatFile;
1143 bM = RmsfFile || bProj;
1144 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1145 || TwoDPlotFile || ThreeDPlotFile;
1146 bIndex = bM || bProj;
1147 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1148 FilterFile || (bIndex && indexfile);
1149 bCompare = Vec2File || Eig2File;
1150 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1152 read_eigenvectors(VecFile, &natoms, &bFit1,
1153 &xref1, &bDMR1, &xav1, &bDMA1,
1154 &nvec1, &eignr1, &eigvec1, &eigval1);
1157 /* Overwrite eigenvalues from separate files if the user provides them */
1158 if (EigFile != NULL)
1160 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1161 if (neig_tmp != neig1)
1163 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1166 srenew(eigval1, neig1);
1167 for (j = 0; j < neig1; j++)
1169 real tmp = eigval1[j];
1170 eigval1[j] = xvgdata[1][j];
1171 if (debug && (eigval1[j] != tmp))
1173 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1174 j, tmp, eigval1[j]);
1177 for (j = 0; j < i; j++)
1182 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1189 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1191 calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1192 calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1199 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1201 read_eigenvectors(Vec2File, &neig2, &bFit2,
1202 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1207 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1211 if (Eig2File != NULL)
1213 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1214 srenew(eigval2, neig2);
1215 for (j = 0; j < neig2; j++)
1217 eigval2[j] = xvgdata[1][j];
1219 for (j = 0; j < i; j++)
1224 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1228 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1232 if ((xref1 == NULL) && (bM || bTraj))
1248 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1249 title, &top, &ePBC, &xtop, NULL, topbox, bM);
1251 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1252 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1253 /* Fitting is only required for the projection */
1258 printf("\nNote: the structure in %s should be the same\n"
1259 " as the one used for the fit in g_covar\n", topfile);
1261 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1262 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1264 snew(w_rls, atoms->nr);
1265 for (i = 0; (i < nfit); i++)
1269 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1273 w_rls[ifit[i]] = 1.0;
1277 snew(xrefp, atoms->nr);
1280 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1283 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);
1285 for (i = 0; (i < nfit); i++)
1287 copy_rvec(xref1[i], xrefp[ifit[i]]);
1292 /* The top coordinates are the fitting reference */
1293 for (i = 0; (i < nfit); i++)
1295 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1297 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1300 gmx_rmpbc_done(gpbc);
1305 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1306 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1309 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1314 snew(sqrtm, natoms);
1317 proj_unit = "u\\S1/2\\Nnm";
1318 for (i = 0; (i < natoms); i++)
1320 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1326 for (i = 0; (i < natoms); i++)
1336 for (i = 0; (i < natoms); i++)
1338 for (d = 0; (d < DIM); d++)
1340 t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1341 totmass += sqr(sqrtm[i]);
1344 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1345 " %.3f (nm)\n\n", sqrt(t/totmass));
1356 /* make an index from first to last */
1357 nout = last-first+1;
1359 for (i = 0; i < nout; i++)
1361 iout[i] = first-1+i;
1364 else if (ThreeDPlotFile)
1366 /* make an index of first+(0,1,2) and last */
1367 nout = bPDB3D ? 4 : 3;
1368 nout = min(last-first+1, nout);
1376 iout[nout-1] = last-1;
1380 /* make an index of first and last */
1389 printf("Select eigenvectors for output, end your selection with 0\n");
1396 srenew(iout, nout+1);
1397 if (1 != scanf("%d", &iout[nout]))
1399 gmx_fatal(FARGS, "Error reading user input");
1403 while (iout[nout] >= 0);
1407 /* make an index of the eigenvectors which are present */
1410 for (i = 0; i < nout; i++)
1413 while ((j < nvec1) && (eignr1[j] != iout[i]))
1417 if ((j < nvec1) && (eignr1[j] == iout[i]))
1419 outvec[noutvec] = j;
1423 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1426 fprintf(stderr, ":");
1427 for (j = 0; j < noutvec; j++)
1429 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1432 fprintf(stderr, "\n");
1436 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1441 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1447 project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1448 bTop ? &top : NULL, ePBC, topbox,
1449 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1450 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1451 bFit1, xrefp, nfit, ifit, w_rls,
1452 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1458 overlap(OverlapFile, natoms,
1459 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1464 inprod_matrix(InpMatFile, natoms,
1465 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1466 bFirstLastSet, noutvec, outvec);
1471 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1475 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1476 !bCompare && !bEntropy)
1478 fprintf(stderr, "\nIf you want some output,"
1479 " set one (or two or ...) of the output file options\n");
1483 view_all(oenv, NFILE, fnm);