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 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
264 fprintf(out, "%10.4f %10.5f\n",
265 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
267 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
274 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
275 real *eigval1, int neig1, real *eigval2, int neig2)
279 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
283 n = min(n, min(neig1, neig2));
284 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
287 for (i = 0; i < n; i++)
294 eigval1[i] = sqrt(eigval1[i]);
297 for (i = n; i < neig1; i++)
299 trace1 += eigval1[i];
302 for (i = 0; i < n; i++)
309 eigval2[i] = sqrt(eigval2[i]);
312 for (i = n; i < neig2; i++)
314 trace2 += eigval2[i];
317 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
318 if (neig1 != n || neig2 != n)
320 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
321 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
323 fprintf(stdout, "Square root of the traces: %g and %g\n",
324 sqrt(sum1), sqrt(sum2));
327 for (i = 0; i < n; i++)
330 for (j = 0; j < n; j++)
333 for (k = 0; k < natoms; k++)
335 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
337 tmp += eigval2[j]*ip*ip;
339 sab += eigval1[i]*tmp;
342 samsb2 = sum1+sum2-2*sab;
348 fprintf(stdout, "The overlap of the covariance matrices:\n");
349 fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
350 tmp = 1-sab/sqrt(sum1*sum2);
355 fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
359 static void inprod_matrix(const char *matfile, int natoms,
360 int nvec1, int *eignr1, rvec **eigvec1,
361 int nvec2, int *eignr2, rvec **eigvec2,
362 gmx_bool bSelect, int noutvec, int *outvec)
366 int i, x1, y1, x, y, nlevels;
368 real inp, *t_x, *t_y, max;
376 for (y1 = 0; y1 < nx; y1++)
378 if (outvec[y1] < nvec2)
380 t_y[ny] = eignr2[outvec[y1]]+1;
389 for (y = 0; y < ny; y++)
391 t_y[y] = eignr2[y]+1;
395 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
401 for (x1 = 0; x1 < nx; x1++)
412 t_x[x1] = eignr1[x]+1;
413 fprintf(stderr, " %d", eignr1[x]+1);
414 for (y1 = 0; y1 < ny; y1++)
419 while (outvec[y1] >= nvec2)
429 for (i = 0; i < natoms; i++)
431 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
433 mat[x1][y1] = fabs(inp);
434 if (mat[x1][y1] > max)
440 fprintf(stderr, "\n");
441 rlo.r = 1; rlo.g = 1; rlo.b = 1;
442 rhi.r = 0; rhi.g = 0; rhi.b = 0;
444 out = gmx_ffopen(matfile, "w");
445 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
446 nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
450 static void overlap(const char *outfile, int natoms,
452 int nvec2, int *eignr2, rvec **eigvec2,
453 int noutvec, int *outvec,
454 const output_env_t oenv)
460 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
461 for (i = 0; i < noutvec; i++)
463 fprintf(stderr, "%d ", outvec[i]+1);
465 fprintf(stderr, "\n");
467 out = xvgropen(outfile, "Subspace overlap",
468 "Eigenvectors of trajectory 2", "Overlap", oenv);
469 if (output_env_get_print_xvgr_codes(oenv))
471 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
474 for (x = 0; x < nvec2; x++)
476 for (v = 0; v < noutvec; v++)
480 for (i = 0; i < natoms; i++)
482 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
486 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
492 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
493 const char *projfile, const char *twodplotfile,
494 const char *threedplotfile, const char *filterfile, int skip,
495 const char *extremefile, gmx_bool bExtrAll, real extreme,
496 int nextr, t_atoms *atoms, int natoms, atom_id *index,
497 gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
498 real *sqrtm, rvec *xav,
499 int *eignr, rvec **eigvec,
500 int noutvec, int *outvec, gmx_bool bSplit,
501 const output_env_t oenv)
503 FILE *xvgrout = NULL;
504 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
505 t_trxstatus *out = NULL;
507 int noutvec_extr, imin, imax;
512 real t, inp, **inprod = NULL, min = 0, max = 0;
513 char str[STRLEN], str2[STRLEN], **ylabel, *c;
515 gmx_rmpbc_t gpbc = NULL;
521 noutvec_extr = noutvec;
531 snew(inprod, noutvec+1);
535 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
537 for (i = 0; i < noutvec; i++)
539 fprintf(stderr, "%d ", outvec[i]+1);
541 fprintf(stderr, "\n");
542 out = open_trx(filterfile, "w");
547 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
550 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);
556 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
559 for (i = 0; i < nat; i++)
569 gmx_rmpbc(gpbc, nat, box, xread);
571 if (nframes >= snew_size)
574 for (i = 0; i < noutvec+1; i++)
576 srenew(inprod[i], snew_size);
579 inprod[noutvec][nframes] = t;
580 /* calculate x: a fitted struture of the selected atoms */
583 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
584 do_fit(nat, w_rls, xref, xread);
586 for (i = 0; i < natoms; i++)
588 copy_rvec(xread[index[i]], x[i]);
591 for (v = 0; v < noutvec; v++)
594 /* calculate (mass-weighted) projection */
596 for (i = 0; i < natoms; i++)
598 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
599 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
600 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
602 inprod[v][nframes] = inp;
606 for (i = 0; i < natoms; i++)
608 for (d = 0; d < DIM; d++)
610 /* misuse xread for output */
611 xread[index[i]][d] = xav[i][d];
612 for (v = 0; v < noutvec; v++)
614 xread[index[i]][d] +=
615 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
619 write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
625 while (read_next_x(oenv, status, &t, xread, box));
635 snew(xread, atoms->nr);
640 gmx_rmpbc_done(gpbc);
646 snew(ylabel, noutvec);
647 for (v = 0; v < noutvec; v++)
649 sprintf(str, "vec %d", eignr[outvec[v]]+1);
650 ylabel[v] = strdup(str);
652 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
653 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
654 (const char **)ylabel,
655 nframes, inprod[noutvec], inprod, NULL,
656 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
661 sprintf(str, "projection on eigenvector %d (%s)",
662 eignr[outvec[0]]+1, proj_unit);
663 sprintf(str2, "projection on eigenvector %d (%s)",
664 eignr[outvec[noutvec-1]]+1, proj_unit);
665 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
667 for (i = 0; i < nframes; i++)
669 if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
671 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
673 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
675 gmx_ffclose(xvgrout);
690 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
694 bPDB = fn2ftp(threedplotfile) == efPDB;
696 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
698 b4D = bPDB && (noutvec >= 4);
701 fprintf(stderr, "You have selected four or more eigenvectors:\n"
702 "fourth eigenvector will be plotted "
703 "in bfactor field of pdb file\n");
704 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
705 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
706 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
710 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
711 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
713 init_t_atoms(&atoms, nframes, FALSE);
717 resnm = strdup("PRJ");
721 fact = 10000.0/nframes;
728 for (i = 0; i < nframes; i++)
730 atoms.atomname[i] = &atnm;
731 atoms.atom[i].resind = i;
732 atoms.resinfo[i].name = &resnm;
733 atoms.resinfo[i].nr = ceil(i*fact);
734 atoms.resinfo[i].ic = ' ';
735 x[i][XX] = inprod[0][i];
736 x[i][YY] = inprod[1][i];
737 x[i][ZZ] = inprod[2][i];
743 if ( ( b4D || bSplit ) && bPDB)
745 out = gmx_ffopen(threedplotfile, "w");
746 fprintf(out, "HEADER %s\n", str);
749 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
752 for (i = 0; i < atoms.nr; i++)
754 if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
756 fprintf(out, "TER\n");
759 gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
760 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
763 fprintf(out, "CONECT%5d%5d\n", i, i+1);
767 fprintf(out, "TER\n");
772 write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
774 free_t_atoms(&atoms, FALSE);
779 snew(pmin, noutvec_extr);
780 snew(pmax, noutvec_extr);
783 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
785 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
788 for (v = 0; v < noutvec_extr; v++)
790 for (i = 0; i < nframes; i++)
792 if (inprod[v][i] < inprod[v][imin])
796 if (inprod[v][i] > inprod[v][imax])
801 pmin[v] = inprod[v][imin];
802 pmax[v] = inprod[v][imax];
803 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
805 pmin[v], imin, pmax[v], imax);
813 /* build format string for filename: */
814 strcpy(str, extremefile); /* copy filename */
815 c = strrchr(str, '.'); /* find where extention begins */
816 strcpy(str2, c); /* get extention */
817 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
818 for (v = 0; v < noutvec_extr; v++)
820 /* make filename using format string */
821 if (noutvec_extr == 1)
823 strcpy(str2, extremefile);
827 sprintf(str2, str, eignr[outvec[v]]+1);
829 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
830 nextr, outvec[v]+1, str2);
831 out = open_trx(str2, "w");
832 for (frame = 0; frame < nextr; frame++)
834 if ((extreme == 0) && (nextr <= 3))
836 for (i = 0; i < natoms; i++)
838 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
841 for (i = 0; i < natoms; i++)
843 for (d = 0; d < DIM; d++)
846 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
847 *eigvec[outvec[v]][i][d]/sqrtm[i]);
850 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
857 fprintf(stderr, "\n");
860 static void components(const char *outfile, int natoms,
861 int *eignr, rvec **eigvec,
862 int noutvec, int *outvec,
863 const output_env_t oenv)
867 char str[STRLEN], **ylabel;
869 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
871 snew(ylabel, noutvec);
874 for (i = 0; i < natoms; i++)
878 for (g = 0; g < noutvec; g++)
881 sprintf(str, "vec %d", eignr[v]+1);
882 ylabel[g] = strdup(str);
884 for (s = 0; s < 4; s++)
886 snew(y[g][s], natoms);
888 for (i = 0; i < natoms; i++)
890 y[g][0][i] = norm(eigvec[v][i]);
891 for (s = 0; s < 3; s++)
893 y[g][s+1][i] = eigvec[v][i][s];
897 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
898 "black: total, red: x, green: y, blue: z",
899 "Atom number", (const char **)ylabel,
900 natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
901 fprintf(stderr, "\n");
904 static void rmsf(const char *outfile, int natoms, real *sqrtm,
905 int *eignr, rvec **eigvec,
906 int noutvec, int *outvec,
907 real *eigval, int neig,
908 const output_env_t oenv)
912 char str[STRLEN], **ylabel;
914 for (i = 0; i < neig; i++)
922 fprintf(stderr, "Writing rmsf to %s\n", outfile);
924 snew(ylabel, noutvec);
927 for (i = 0; i < natoms; i++)
931 for (g = 0; g < noutvec; g++)
934 if (eignr[v] >= neig)
936 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
938 sprintf(str, "vec %d", eignr[v]+1);
939 ylabel[g] = strdup(str);
941 for (i = 0; i < natoms; i++)
943 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
946 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
947 "Atom number", (const char **)ylabel,
948 natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
949 fprintf(stderr, "\n");
952 int gmx_anaeig(int argc, char *argv[])
954 static const char *desc[] = {
955 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
956 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
957 "([gmx-nmeig]).[PAR]",
959 "When a trajectory is projected on eigenvectors, all structures are",
960 "fitted to the structure in the eigenvector file, if present, otherwise",
961 "to the structure in the structure file. When no run input file is",
962 "supplied, periodicity will not be taken into account. Most analyses",
963 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
964 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
966 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
967 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
969 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
970 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
972 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
973 "[TT]-first[tt] to [TT]-last[tt].",
974 "The projections of a trajectory on the eigenvectors of its",
975 "covariance matrix are called principal components (pc's).",
976 "It is often useful to check the cosine content of the pc's,",
977 "since the pc's of random diffusion are cosines with the number",
978 "of periods equal to half the pc index.",
979 "The cosine content of the pc's can be calculated with the program",
980 "[gmx-analyze].[PAR]",
982 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
983 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
985 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
986 "three selected eigenvectors.[PAR]",
988 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
989 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
991 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
992 "on the average structure and interpolate [TT]-nframes[tt] frames",
993 "between them, or set your own extremes with [TT]-max[tt]. The",
994 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
995 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
996 "will be written to separate files. Chain identifiers will be added",
997 "when writing a [TT].pdb[tt] file with two or three structures (you",
998 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1000 " Overlap calculations between covariance analysis:[BR]",
1001 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1003 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1004 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1005 "in file [TT]-v[tt].[PAR]",
1007 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1008 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1009 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1010 "have been set explicitly.[PAR]",
1012 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1013 "a single number for the overlap between the covariance matrices is",
1014 "generated. The formulas are:[BR]",
1015 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1016 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1017 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1018 "where M1 and M2 are the two covariance matrices and tr is the trace",
1019 "of a matrix. The numbers are proportional to the overlap of the square",
1020 "root of the fluctuations. The normalized overlap is the most useful",
1021 "number, it is 1 for identical matrices and 0 when the sampled",
1022 "subspaces are orthogonal.[PAR]",
1023 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1024 "computed based on the Quasiharmonic approach and based on",
1025 "Schlitter's formula."
1027 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1028 static real max = 0.0, temp = 298.15;
1029 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1031 { "-first", FALSE, etINT, {&first},
1032 "First eigenvector for analysis (-1 is select)" },
1033 { "-last", FALSE, etINT, {&last},
1034 "Last eigenvector for analysis (-1 is till the last)" },
1035 { "-skip", FALSE, etINT, {&skip},
1036 "Only analyse every nr-th frame" },
1037 { "-max", FALSE, etREAL, {&max},
1038 "Maximum for projection of the eigenvector on the average structure, "
1039 "max=0 gives the extremes" },
1040 { "-nframes", FALSE, etINT, {&nextr},
1041 "Number of frames for the extremes output" },
1042 { "-split", FALSE, etBOOL, {&bSplit},
1043 "Split eigenvector projections where time is zero" },
1044 { "-entropy", FALSE, etBOOL, {&bEntropy},
1045 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1046 { "-temp", FALSE, etREAL, {&temp},
1047 "Temperature for entropy calculations" },
1048 { "-nevskip", FALSE, etINT, {&nskip},
1049 "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." }
1051 #define NPA asize(pa)
1057 t_atoms *atoms = NULL;
1058 rvec *xtop, *xref1, *xref2, *xrefp = NULL;
1059 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1060 int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1061 rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1063 real xid, totmass, *sqrtm, *w_rls, t, lambda;
1066 const char *indexfile;
1069 int nout, *iout, noutvec, *outvec, nfit;
1070 atom_id *index, *ifit;
1071 const char *VecFile, *Vec2File, *topfile;
1072 const char *EigFile, *Eig2File;
1073 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1074 const char *TwoDPlotFile, *ThreeDPlotFile;
1075 const char *FilterFile, *ExtremeFile;
1076 const char *OverlapFile, *InpMatFile;
1077 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1078 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1079 real *eigval1 = NULL, *eigval2 = NULL;
1086 { efTRN, "-v", "eigenvec", ffREAD },
1087 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1088 { efTRX, "-f", NULL, ffOPTRD },
1089 { efTPS, NULL, NULL, ffOPTRD },
1090 { efNDX, NULL, NULL, ffOPTRD },
1091 { efXVG, "-eig", "eigenval", ffOPTRD },
1092 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1093 { efXVG, "-comp", "eigcomp", ffOPTWR },
1094 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1095 { efXVG, "-proj", "proj", ffOPTWR },
1096 { efXVG, "-2d", "2dproj", ffOPTWR },
1097 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1098 { efTRX, "-filt", "filtered", ffOPTWR },
1099 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1100 { efXVG, "-over", "overlap", ffOPTWR },
1101 { efXPM, "-inpr", "inprod", ffOPTWR }
1103 #define NFILE asize(fnm)
1105 if (!parse_common_args(&argc, argv,
1106 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
1107 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1112 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1114 VecFile = opt2fn("-v", NFILE, fnm);
1115 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1116 topfile = ftp2fn(efTPS, NFILE, fnm);
1117 EigFile = opt2fn_null("-eig", NFILE, fnm);
1118 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1119 CompFile = opt2fn_null("-comp", NFILE, fnm);
1120 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1121 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1122 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1123 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1124 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1125 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1126 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1127 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1129 bTop = fn2bTPX(topfile);
1130 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1131 || FilterFile || ExtremeFile;
1133 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1134 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1135 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1136 bVec2 = Vec2File || OverlapFile || InpMatFile;
1137 bM = RmsfFile || bProj;
1138 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1139 || TwoDPlotFile || ThreeDPlotFile;
1140 bIndex = bM || bProj;
1141 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1142 FilterFile || (bIndex && indexfile);
1143 bCompare = Vec2File || Eig2File;
1144 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1146 read_eigenvectors(VecFile, &natoms, &bFit1,
1147 &xref1, &bDMR1, &xav1, &bDMA1,
1148 &nvec1, &eignr1, &eigvec1, &eigval1);
1151 /* Overwrite eigenvalues from separate files if the user provides them */
1152 if (EigFile != NULL)
1154 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1155 if (neig_tmp != neig1)
1157 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1160 srenew(eigval1, neig1);
1161 for (j = 0; j < neig1; j++)
1163 real tmp = eigval1[j];
1164 eigval1[j] = xvgdata[1][j];
1165 if (debug && (eigval1[j] != tmp))
1167 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1168 j, tmp, eigval1[j]);
1171 for (j = 0; j < i; j++)
1176 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1183 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1185 calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1186 calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1193 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1195 read_eigenvectors(Vec2File, &neig2, &bFit2,
1196 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1201 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1205 if (Eig2File != NULL)
1207 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1208 srenew(eigval2, neig2);
1209 for (j = 0; j < neig2; j++)
1211 eigval2[j] = xvgdata[1][j];
1213 for (j = 0; j < i; j++)
1218 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1222 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1226 if ((xref1 == NULL) && (bM || bTraj))
1242 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1243 title, &top, &ePBC, &xtop, NULL, topbox, bM);
1245 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1246 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1247 /* Fitting is only required for the projection */
1252 printf("\nNote: the structure in %s should be the same\n"
1253 " as the one used for the fit in g_covar\n", topfile);
1255 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1256 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1258 snew(w_rls, atoms->nr);
1259 for (i = 0; (i < nfit); i++)
1263 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1267 w_rls[ifit[i]] = 1.0;
1271 snew(xrefp, atoms->nr);
1274 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1277 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);
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);