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"
45 #include "gromacs/utility/smalloc.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"
65 #include "gromacs/math/do_fit.h"
67 static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
70 double hwkT, w, dS, S = 0;
73 hbar = PLANCK1/(2*M_PI);
74 for (i = 0; (i < n-nskip); i++)
78 lambda = eigval[i]*AMU;
79 w = sqrt(BOLTZMANN*temp/lambda)/NANO;
80 hwkT = (hbar*w)/(BOLTZMANN*temp);
81 dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
85 fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
86 i, w, lambda, hwkT, dS);
91 fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
95 fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
99 static void calc_entropy_schlitter(FILE *fp, int n, int nskip,
100 real *eigval, real temp)
106 double hbar, kt, kteh, S;
108 hbar = PLANCK1/(2*M_PI);
110 kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
113 fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
117 for (i = 0; (i < n-nskip); i++)
119 dd = 1+kteh*eigval[i];
124 fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S);
127 const char *proj_unit;
129 static real tick_spacing(real range, int minticks)
138 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
139 while (range/sp < minticks-1)
147 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
148 const char *title, const char *subtitle,
149 const char *xlabel, const char **ylabel,
150 int n, real *x, real **y, real ***sy,
151 real scale_x, gmx_bool bZero, gmx_bool bSplit,
152 const output_env_t oenv)
156 real min, max, xsp, ysp;
158 out = gmx_ffopen(file, "w");
159 if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
161 fprintf(out, "@ autoscale onread none\n");
163 for (g = 0; g < ngraphs; g++)
169 for (i = 0; i < n; i++)
185 for (s = 0; s < nsetspergraph; s++)
187 for (i = 0; i < n; i++)
189 if (sy[g][s][i] < min)
193 if (sy[g][s][i] > max)
206 min = min-0.1*(max-min);
208 max = max+0.1*(max-min);
209 xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
210 ysp = tick_spacing(max-min, 3);
211 if (output_env_get_print_xvgr_codes(oenv))
213 fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
216 fprintf(out, "@ title \"%s\"\n", title);
219 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
224 fprintf(out, "@ xaxis label \"%s\"\n", xlabel);
228 fprintf(out, "@ xaxis ticklabel off\n");
232 fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
233 fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
234 fprintf(out, "@ world ymin %g\n", min);
235 fprintf(out, "@ world ymax %g\n", max);
237 fprintf(out, "@ view xmin 0.15\n");
238 fprintf(out, "@ view xmax 0.85\n");
239 fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs);
240 fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs);
241 fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
242 fprintf(out, "@ xaxis tick major %g\n", xsp);
243 fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
244 fprintf(out, "@ xaxis ticklabel start type spec\n");
245 fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp);
246 fprintf(out, "@ yaxis tick major %g\n", ysp);
247 fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
248 fprintf(out, "@ yaxis ticklabel start type spec\n");
249 fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp);
250 if ((min < 0) && (max > 0))
252 fprintf(out, "@ zeroxaxis bar on\n");
253 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
256 for (s = 0; s < nsetspergraph; s++)
258 for (i = 0; i < n; i++)
260 if (bSplit && i > 0 && fabs(x[i]) < 1e-5)
262 if (output_env_get_print_xvgr_codes(oenv))
267 fprintf(out, "%10.4f %10.5f\n",
268 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
270 if (output_env_get_print_xvgr_codes(oenv))
280 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
281 real *eigval1, int neig1, real *eigval2, int neig2)
285 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
289 n = min(n, min(neig1, neig2));
290 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
293 for (i = 0; i < n; i++)
300 eigval1[i] = sqrt(eigval1[i]);
303 for (i = n; i < neig1; i++)
305 trace1 += eigval1[i];
308 for (i = 0; i < n; i++)
315 eigval2[i] = sqrt(eigval2[i]);
318 for (i = n; i < neig2; i++)
320 trace2 += eigval2[i];
323 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
324 if (neig1 != n || neig2 != n)
326 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
327 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
329 fprintf(stdout, "Square root of the traces: %g and %g\n",
330 sqrt(sum1), sqrt(sum2));
333 for (i = 0; i < n; i++)
336 for (j = 0; j < n; j++)
339 for (k = 0; k < natoms; k++)
341 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
343 tmp += eigval2[j]*ip*ip;
345 sab += eigval1[i]*tmp;
348 samsb2 = sum1+sum2-2*sab;
354 fprintf(stdout, "The overlap of the covariance matrices:\n");
355 fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
356 tmp = 1-sab/sqrt(sum1*sum2);
361 fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
365 static void inprod_matrix(const char *matfile, int natoms,
366 int nvec1, int *eignr1, rvec **eigvec1,
367 int nvec2, int *eignr2, rvec **eigvec2,
368 gmx_bool bSelect, int noutvec, int *outvec)
372 int i, x1, y1, x, y, nlevels;
374 real inp, *t_x, *t_y, max;
382 for (y1 = 0; y1 < nx; y1++)
384 if (outvec[y1] < nvec2)
386 t_y[ny] = eignr2[outvec[y1]]+1;
395 for (y = 0; y < ny; y++)
397 t_y[y] = eignr2[y]+1;
401 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
407 for (x1 = 0; x1 < nx; x1++)
418 t_x[x1] = eignr1[x]+1;
419 fprintf(stderr, " %d", eignr1[x]+1);
420 for (y1 = 0; y1 < ny; y1++)
425 while (outvec[y1] >= nvec2)
435 for (i = 0; i < natoms; i++)
437 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
439 mat[x1][y1] = fabs(inp);
440 if (mat[x1][y1] > max)
446 fprintf(stderr, "\n");
447 rlo.r = 1; rlo.g = 1; rlo.b = 1;
448 rhi.r = 0; rhi.g = 0; rhi.b = 0;
450 out = gmx_ffopen(matfile, "w");
451 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
452 nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
456 static void overlap(const char *outfile, int natoms,
458 int nvec2, int *eignr2, rvec **eigvec2,
459 int noutvec, int *outvec,
460 const output_env_t oenv)
466 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
467 for (i = 0; i < noutvec; i++)
469 fprintf(stderr, "%d ", outvec[i]+1);
471 fprintf(stderr, "\n");
473 out = xvgropen(outfile, "Subspace overlap",
474 "Eigenvectors of trajectory 2", "Overlap", oenv);
475 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n",
478 for (x = 0; x < nvec2; x++)
480 for (v = 0; v < noutvec; v++)
484 for (i = 0; i < natoms; i++)
486 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
490 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
496 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
497 const char *projfile, const char *twodplotfile,
498 const char *threedplotfile, const char *filterfile, int skip,
499 const char *extremefile, gmx_bool bExtrAll, real extreme,
500 int nextr, t_atoms *atoms, int natoms, atom_id *index,
501 gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
502 real *sqrtm, rvec *xav,
503 int *eignr, rvec **eigvec,
504 int noutvec, int *outvec, gmx_bool bSplit,
505 const output_env_t oenv)
507 FILE *xvgrout = NULL;
508 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
509 t_trxstatus *out = NULL;
511 int noutvec_extr, imin, imax;
516 real t, inp, **inprod = NULL, min = 0, max = 0;
517 char str[STRLEN], str2[STRLEN], **ylabel, *c;
519 gmx_rmpbc_t gpbc = NULL;
525 noutvec_extr = noutvec;
535 snew(inprod, noutvec+1);
539 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
541 for (i = 0; i < noutvec; i++)
543 fprintf(stderr, "%d ", outvec[i]+1);
545 fprintf(stderr, "\n");
546 out = open_trx(filterfile, "w");
551 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
554 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);
560 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
563 for (i = 0; i < nat; i++)
573 gmx_rmpbc(gpbc, nat, box, xread);
575 if (nframes >= snew_size)
578 for (i = 0; i < noutvec+1; i++)
580 srenew(inprod[i], snew_size);
583 inprod[noutvec][nframes] = t;
584 /* calculate x: a fitted struture of the selected atoms */
587 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
588 do_fit(nat, w_rls, xref, xread);
590 for (i = 0; i < natoms; i++)
592 copy_rvec(xread[index[i]], x[i]);
595 for (v = 0; v < noutvec; v++)
598 /* calculate (mass-weighted) projection */
600 for (i = 0; i < natoms; i++)
602 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
603 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
604 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
606 inprod[v][nframes] = inp;
610 for (i = 0; i < natoms; i++)
612 for (d = 0; d < DIM; d++)
614 /* misuse xread for output */
615 xread[index[i]][d] = xav[i][d];
616 for (v = 0; v < noutvec; v++)
618 xread[index[i]][d] +=
619 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
623 write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
629 while (read_next_x(oenv, status, &t, xread, box));
639 snew(xread, atoms->nr);
644 gmx_rmpbc_done(gpbc);
650 snew(ylabel, noutvec);
651 for (v = 0; v < noutvec; v++)
653 sprintf(str, "vec %d", eignr[outvec[v]]+1);
654 ylabel[v] = strdup(str);
656 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
657 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
658 (const char **)ylabel,
659 nframes, inprod[noutvec], inprod, NULL,
660 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
665 sprintf(str, "projection on eigenvector %d (%s)",
666 eignr[outvec[0]]+1, proj_unit);
667 sprintf(str2, "projection on eigenvector %d (%s)",
668 eignr[outvec[noutvec-1]]+1, proj_unit);
669 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
671 for (i = 0; i < nframes; i++)
673 if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
675 fprintf(xvgrout, "&\n");
677 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
679 gmx_ffclose(xvgrout);
688 char *resnm, *atnm, pdbform[STRLEN];
694 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
698 bPDB = fn2ftp(threedplotfile) == efPDB;
700 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
702 b4D = bPDB && (noutvec >= 4);
705 fprintf(stderr, "You have selected four or more eigenvectors:\n"
706 "fourth eigenvector will be plotted "
707 "in bfactor field of pdb file\n");
708 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
709 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
710 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
714 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
715 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
717 init_t_atoms(&atoms, nframes, FALSE);
721 resnm = strdup("PRJ");
725 fact = 10000.0/nframes;
732 for (i = 0; i < nframes; i++)
734 atoms.atomname[i] = &atnm;
735 atoms.atom[i].resind = i;
736 atoms.resinfo[i].name = &resnm;
737 atoms.resinfo[i].nr = ceil(i*fact);
738 atoms.resinfo[i].ic = ' ';
739 x[i][XX] = inprod[0][i];
740 x[i][YY] = inprod[1][i];
741 x[i][ZZ] = inprod[2][i];
747 if ( ( b4D || bSplit ) && bPDB)
749 strcpy(pdbform, get_pdbformat());
750 strcat(pdbform, "%8.4f%8.4f\n");
752 out = gmx_ffopen(threedplotfile, "w");
753 fprintf(out, "HEADER %s\n", str);
756 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
759 for (i = 0; i < atoms.nr; i++)
761 if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
763 fprintf(out, "TER\n");
766 fprintf(out, pdbform, "ATOM", i+1, "C", "PRJ", ' ', j+1,
767 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i]);
770 fprintf(out, "CONECT%5d%5d\n", i, i+1);
774 fprintf(out, "TER\n");
779 write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
781 free_t_atoms(&atoms, FALSE);
786 snew(pmin, noutvec_extr);
787 snew(pmax, noutvec_extr);
790 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
792 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
795 for (v = 0; v < noutvec_extr; v++)
797 for (i = 0; i < nframes; i++)
799 if (inprod[v][i] < inprod[v][imin])
803 if (inprod[v][i] > inprod[v][imax])
808 pmin[v] = inprod[v][imin];
809 pmax[v] = inprod[v][imax];
810 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
812 pmin[v], imin, pmax[v], imax);
820 /* build format string for filename: */
821 strcpy(str, extremefile); /* copy filename */
822 c = strrchr(str, '.'); /* find where extention begins */
823 strcpy(str2, c); /* get extention */
824 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
825 for (v = 0; v < noutvec_extr; v++)
827 /* make filename using format string */
828 if (noutvec_extr == 1)
830 strcpy(str2, extremefile);
834 sprintf(str2, str, eignr[outvec[v]]+1);
836 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
837 nextr, outvec[v]+1, str2);
838 out = open_trx(str2, "w");
839 for (frame = 0; frame < nextr; frame++)
841 if ((extreme == 0) && (nextr <= 3))
843 for (i = 0; i < natoms; i++)
845 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
848 for (i = 0; i < natoms; i++)
850 for (d = 0; d < DIM; d++)
853 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
854 *eigvec[outvec[v]][i][d]/sqrtm[i]);
857 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
864 fprintf(stderr, "\n");
867 static void components(const char *outfile, int natoms,
868 int *eignr, rvec **eigvec,
869 int noutvec, int *outvec,
870 const output_env_t oenv)
874 char str[STRLEN], **ylabel;
876 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
878 snew(ylabel, noutvec);
881 for (i = 0; i < natoms; i++)
885 for (g = 0; g < noutvec; g++)
888 sprintf(str, "vec %d", eignr[v]+1);
889 ylabel[g] = strdup(str);
891 for (s = 0; s < 4; s++)
893 snew(y[g][s], natoms);
895 for (i = 0; i < natoms; i++)
897 y[g][0][i] = norm(eigvec[v][i]);
898 for (s = 0; s < 3; s++)
900 y[g][s+1][i] = eigvec[v][i][s];
904 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
905 "black: total, red: x, green: y, blue: z",
906 "Atom number", (const char **)ylabel,
907 natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
908 fprintf(stderr, "\n");
911 static void rmsf(const char *outfile, int natoms, real *sqrtm,
912 int *eignr, rvec **eigvec,
913 int noutvec, int *outvec,
914 real *eigval, int neig,
915 const output_env_t oenv)
919 char str[STRLEN], **ylabel;
921 for (i = 0; i < neig; i++)
929 fprintf(stderr, "Writing rmsf to %s\n", outfile);
931 snew(ylabel, noutvec);
934 for (i = 0; i < natoms; i++)
938 for (g = 0; g < noutvec; g++)
941 if (eignr[v] >= neig)
943 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
945 sprintf(str, "vec %d", eignr[v]+1);
946 ylabel[g] = strdup(str);
948 for (i = 0; i < natoms; i++)
950 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
953 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
954 "Atom number", (const char **)ylabel,
955 natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
956 fprintf(stderr, "\n");
959 int gmx_anaeig(int argc, char *argv[])
961 static const char *desc[] = {
962 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
963 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
964 "([gmx-nmeig]).[PAR]",
966 "When a trajectory is projected on eigenvectors, all structures are",
967 "fitted to the structure in the eigenvector file, if present, otherwise",
968 "to the structure in the structure file. When no run input file is",
969 "supplied, periodicity will not be taken into account. Most analyses",
970 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
971 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
973 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
974 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
976 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
977 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
979 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
980 "[TT]-first[tt] to [TT]-last[tt].",
981 "The projections of a trajectory on the eigenvectors of its",
982 "covariance matrix are called principal components (pc's).",
983 "It is often useful to check the cosine content of the pc's,",
984 "since the pc's of random diffusion are cosines with the number",
985 "of periods equal to half the pc index.",
986 "The cosine content of the pc's can be calculated with the program",
987 "[gmx-analyze].[PAR]",
989 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
990 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
992 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
993 "three selected eigenvectors.[PAR]",
995 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
996 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
998 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
999 "on the average structure and interpolate [TT]-nframes[tt] frames",
1000 "between them, or set your own extremes with [TT]-max[tt]. The",
1001 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
1002 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1003 "will be written to separate files. Chain identifiers will be added",
1004 "when writing a [TT].pdb[tt] file with two or three structures (you",
1005 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1007 " Overlap calculations between covariance analysis:[BR]",
1008 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1010 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1011 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1012 "in file [TT]-v[tt].[PAR]",
1014 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1015 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1016 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1017 "have been set explicitly.[PAR]",
1019 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1020 "a single number for the overlap between the covariance matrices is",
1021 "generated. The formulas are:[BR]",
1022 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1023 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1024 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1025 "where M1 and M2 are the two covariance matrices and tr is the trace",
1026 "of a matrix. The numbers are proportional to the overlap of the square",
1027 "root of the fluctuations. The normalized overlap is the most useful",
1028 "number, it is 1 for identical matrices and 0 when the sampled",
1029 "subspaces are orthogonal.[PAR]",
1030 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1031 "computed based on the Quasiharmonic approach and based on",
1032 "Schlitter's formula."
1034 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1035 static real max = 0.0, temp = 298.15;
1036 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1038 { "-first", FALSE, etINT, {&first},
1039 "First eigenvector for analysis (-1 is select)" },
1040 { "-last", FALSE, etINT, {&last},
1041 "Last eigenvector for analysis (-1 is till the last)" },
1042 { "-skip", FALSE, etINT, {&skip},
1043 "Only analyse every nr-th frame" },
1044 { "-max", FALSE, etREAL, {&max},
1045 "Maximum for projection of the eigenvector on the average structure, "
1046 "max=0 gives the extremes" },
1047 { "-nframes", FALSE, etINT, {&nextr},
1048 "Number of frames for the extremes output" },
1049 { "-split", FALSE, etBOOL, {&bSplit},
1050 "Split eigenvector projections where time is zero" },
1051 { "-entropy", FALSE, etBOOL, {&bEntropy},
1052 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1053 { "-temp", FALSE, etREAL, {&temp},
1054 "Temperature for entropy calculations" },
1055 { "-nevskip", FALSE, etINT, {&nskip},
1056 "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." }
1058 #define NPA asize(pa)
1064 t_atoms *atoms = NULL;
1065 rvec *xtop, *xref1, *xref2, *xrefp = NULL;
1066 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1067 int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1068 rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1070 real xid, totmass, *sqrtm, *w_rls, t, lambda;
1073 const char *indexfile;
1076 int nout, *iout, noutvec, *outvec, nfit;
1077 atom_id *index, *ifit;
1078 const char *VecFile, *Vec2File, *topfile;
1079 const char *EigFile, *Eig2File;
1080 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1081 const char *TwoDPlotFile, *ThreeDPlotFile;
1082 const char *FilterFile, *ExtremeFile;
1083 const char *OverlapFile, *InpMatFile;
1084 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1085 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1086 real *eigval1 = NULL, *eigval2 = NULL;
1093 { efTRN, "-v", "eigenvec", ffREAD },
1094 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1095 { efTRX, "-f", NULL, ffOPTRD },
1096 { efTPS, NULL, NULL, ffOPTRD },
1097 { efNDX, NULL, NULL, ffOPTRD },
1098 { efXVG, "-eig", "eigenval", ffOPTRD },
1099 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1100 { efXVG, "-comp", "eigcomp", ffOPTWR },
1101 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1102 { efXVG, "-proj", "proj", ffOPTWR },
1103 { efXVG, "-2d", "2dproj", ffOPTWR },
1104 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1105 { efTRX, "-filt", "filtered", ffOPTWR },
1106 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1107 { efXVG, "-over", "overlap", ffOPTWR },
1108 { efXPM, "-inpr", "inprod", ffOPTWR }
1110 #define NFILE asize(fnm)
1112 if (!parse_common_args(&argc, argv,
1113 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
1114 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1119 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1121 VecFile = opt2fn("-v", NFILE, fnm);
1122 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1123 topfile = ftp2fn(efTPS, NFILE, fnm);
1124 EigFile = opt2fn_null("-eig", NFILE, fnm);
1125 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1126 CompFile = opt2fn_null("-comp", NFILE, fnm);
1127 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1128 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1129 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1130 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1131 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1132 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1133 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1134 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1136 bTop = fn2bTPX(topfile);
1137 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1138 || FilterFile || ExtremeFile;
1140 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1141 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1142 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1143 bVec2 = Vec2File || OverlapFile || InpMatFile;
1144 bM = RmsfFile || bProj;
1145 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1146 || TwoDPlotFile || ThreeDPlotFile;
1147 bIndex = bM || bProj;
1148 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1149 FilterFile || (bIndex && indexfile);
1150 bCompare = Vec2File || Eig2File;
1151 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1153 read_eigenvectors(VecFile, &natoms, &bFit1,
1154 &xref1, &bDMR1, &xav1, &bDMA1,
1155 &nvec1, &eignr1, &eigvec1, &eigval1);
1158 /* Overwrite eigenvalues from separate files if the user provides them */
1159 if (EigFile != NULL)
1161 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1162 if (neig_tmp != neig1)
1164 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1167 srenew(eigval1, neig1);
1168 for (j = 0; j < neig1; j++)
1170 real tmp = eigval1[j];
1171 eigval1[j] = xvgdata[1][j];
1172 if (debug && (eigval1[j] != tmp))
1174 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1175 j, tmp, eigval1[j]);
1178 for (j = 0; j < i; j++)
1183 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1190 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1192 calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1193 calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1200 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1202 read_eigenvectors(Vec2File, &neig2, &bFit2,
1203 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1208 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1212 if (Eig2File != NULL)
1214 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1215 srenew(eigval2, neig2);
1216 for (j = 0; j < neig2; j++)
1218 eigval2[j] = xvgdata[1][j];
1220 for (j = 0; j < i; j++)
1225 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1229 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1233 if ((xref1 == NULL) && (bM || bTraj))
1249 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1250 title, &top, &ePBC, &xtop, NULL, topbox, bM);
1252 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1253 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1254 /* Fitting is only required for the projection */
1259 printf("\nNote: the structure in %s should be the same\n"
1260 " as the one used for the fit in g_covar\n", topfile);
1262 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1263 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1265 snew(w_rls, atoms->nr);
1266 for (i = 0; (i < nfit); i++)
1270 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1274 w_rls[ifit[i]] = 1.0;
1278 snew(xrefp, atoms->nr);
1281 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1284 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);
1286 for (i = 0; (i < nfit); i++)
1288 copy_rvec(xref1[i], xrefp[ifit[i]]);
1293 /* The top coordinates are the fitting reference */
1294 for (i = 0; (i < nfit); i++)
1296 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1298 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1301 gmx_rmpbc_done(gpbc);
1306 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1307 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1310 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1315 snew(sqrtm, natoms);
1318 proj_unit = "u\\S1/2\\Nnm";
1319 for (i = 0; (i < natoms); i++)
1321 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1327 for (i = 0; (i < natoms); i++)
1337 for (i = 0; (i < natoms); i++)
1339 for (d = 0; (d < DIM); d++)
1341 t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1342 totmass += sqr(sqrtm[i]);
1345 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1346 " %.3f (nm)\n\n", sqrt(t/totmass));
1357 /* make an index from first to last */
1358 nout = last-first+1;
1360 for (i = 0; i < nout; i++)
1362 iout[i] = first-1+i;
1365 else if (ThreeDPlotFile)
1367 /* make an index of first+(0,1,2) and last */
1368 nout = bPDB3D ? 4 : 3;
1369 nout = min(last-first+1, nout);
1377 iout[nout-1] = last-1;
1381 /* make an index of first and last */
1390 printf("Select eigenvectors for output, end your selection with 0\n");
1397 srenew(iout, nout+1);
1398 if (1 != scanf("%d", &iout[nout]))
1400 gmx_fatal(FARGS, "Error reading user input");
1404 while (iout[nout] >= 0);
1408 /* make an index of the eigenvectors which are present */
1411 for (i = 0; i < nout; i++)
1414 while ((j < nvec1) && (eignr1[j] != iout[i]))
1418 if ((j < nvec1) && (eignr1[j] == iout[i]))
1420 outvec[noutvec] = j;
1424 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1427 fprintf(stderr, ":");
1428 for (j = 0; j < noutvec; j++)
1430 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1433 fprintf(stderr, "\n");
1437 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1442 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1448 project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1449 bTop ? &top : NULL, ePBC, topbox,
1450 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1451 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1452 bFit1, xrefp, nfit, ifit, w_rls,
1453 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1459 overlap(OverlapFile, natoms,
1460 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1465 inprod_matrix(InpMatFile, natoms,
1466 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1467 bFirstLastSet, noutvec, outvec);
1472 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1476 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1477 !bCompare && !bEntropy)
1479 fprintf(stderr, "\nIf you want some output,"
1480 " set one (or two or ...) of the output file options\n");
1484 view_all(oenv, NFILE, fnm);