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.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/legacyheaders/typedefs.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/legacyheaders/macros.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/utility/futil.h"
52 #include "gromacs/topology/index.h"
53 #include "gromacs/fileio/pdbio.h"
54 #include "gromacs/fileio/confio.h"
55 #include "gromacs/fileio/tpxio.h"
56 #include "gromacs/fileio/trxio.h"
57 #include "gromacs/fileio/matio.h"
58 #include "gromacs/fileio/xvgr.h"
59 #include "gromacs/legacyheaders/viewit.h"
60 #include "gromacs/pbcutil/rmpbc.h"
61 #include "gromacs/legacyheaders/txtdump.h"
63 #include "gromacs/math/units.h"
66 #include "gromacs/math/do_fit.h"
68 static void calc_entropy_qh(FILE *fp, int n, real eigval[], real temp, int nskip)
71 double hwkT, w, dS, S = 0;
74 hbar = PLANCK1/(2*M_PI);
75 for (i = 0; (i < n-nskip); i++)
79 lambda = eigval[i]*AMU;
80 w = sqrt(BOLTZMANN*temp/lambda)/NANO;
81 hwkT = (hbar*w)/(BOLTZMANN*temp);
82 dS = (hwkT/(exp(hwkT) - 1) - log(1-exp(-hwkT)));
86 fprintf(debug, "i = %5d w = %10g lam = %10g hwkT = %10g dS = %10g\n",
87 i, w, lambda, hwkT, dS);
92 fprintf(stderr, "eigval[%d] = %g\n", i, eigval[i]);
96 fprintf(fp, "The Entropy due to the Quasi Harmonic approximation is %g J/mol K\n",
100 static void calc_entropy_schlitter(FILE *fp, int n, int nskip,
101 real *eigval, real temp)
107 double hbar, kt, kteh, S;
109 hbar = PLANCK1/(2*M_PI);
111 kteh = kt*exp(2.0)/(hbar*hbar)*AMU*(NANO*NANO);
114 fprintf(debug, "n = %d, nskip = %d kteh = %g\n", n, nskip, kteh);
118 for (i = 0; (i < n-nskip); i++)
120 dd = 1+kteh*eigval[i];
125 fprintf(fp, "The Entropy due to the Schlitter formula is %g J/mol K\n", S);
128 const char *proj_unit;
130 static real tick_spacing(real range, int minticks)
139 sp = 0.2*exp(log(10)*ceil(log(range)/log(10)));
140 while (range/sp < minticks-1)
148 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
149 const char *title, const char *subtitle,
150 const char *xlabel, const char **ylabel,
151 int n, real *x, real **y, real ***sy,
152 real scale_x, gmx_bool bZero, gmx_bool bSplit,
153 const output_env_t oenv)
157 real min, max, xsp, ysp;
159 out = gmx_ffopen(file, "w");
160 if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
162 fprintf(out, "@ autoscale onread none\n");
164 for (g = 0; g < ngraphs; g++)
170 for (i = 0; i < n; i++)
186 for (s = 0; s < nsetspergraph; s++)
188 for (i = 0; i < n; i++)
190 if (sy[g][s][i] < min)
194 if (sy[g][s][i] > max)
207 min = min-0.1*(max-min);
209 max = max+0.1*(max-min);
210 xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
211 ysp = tick_spacing(max-min, 3);
212 if (output_env_get_print_xvgr_codes(oenv))
214 fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
217 fprintf(out, "@ title \"%s\"\n", title);
220 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
225 fprintf(out, "@ xaxis label \"%s\"\n", xlabel);
229 fprintf(out, "@ xaxis ticklabel off\n");
233 fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
234 fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
235 fprintf(out, "@ world ymin %g\n", min);
236 fprintf(out, "@ world ymax %g\n", max);
238 fprintf(out, "@ view xmin 0.15\n");
239 fprintf(out, "@ view xmax 0.85\n");
240 fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs);
241 fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs);
242 fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
243 fprintf(out, "@ xaxis tick major %g\n", xsp);
244 fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
245 fprintf(out, "@ xaxis ticklabel start type spec\n");
246 fprintf(out, "@ xaxis ticklabel start %g\n", ceil(min/xsp)*xsp);
247 fprintf(out, "@ yaxis tick major %g\n", ysp);
248 fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
249 fprintf(out, "@ yaxis ticklabel start type spec\n");
250 fprintf(out, "@ yaxis ticklabel start %g\n", ceil(min/ysp)*ysp);
251 if ((min < 0) && (max > 0))
253 fprintf(out, "@ zeroxaxis bar on\n");
254 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
257 for (s = 0; s < nsetspergraph; s++)
259 for (i = 0; i < n; i++)
261 if (bSplit && i > 0 && fabs(x[i]) < 1e-5)
263 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
265 fprintf(out, "%10.4f %10.5f\n",
266 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
268 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
275 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
276 real *eigval1, int neig1, real *eigval2, int neig2)
280 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
284 n = min(n, min(neig1, neig2));
285 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
288 for (i = 0; i < n; i++)
295 eigval1[i] = sqrt(eigval1[i]);
298 for (i = n; i < neig1; i++)
300 trace1 += eigval1[i];
303 for (i = 0; i < n; i++)
310 eigval2[i] = sqrt(eigval2[i]);
313 for (i = n; i < neig2; i++)
315 trace2 += eigval2[i];
318 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
319 if (neig1 != n || neig2 != n)
321 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
322 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
324 fprintf(stdout, "Square root of the traces: %g and %g\n",
325 sqrt(sum1), sqrt(sum2));
328 for (i = 0; i < n; i++)
331 for (j = 0; j < n; j++)
334 for (k = 0; k < natoms; k++)
336 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
338 tmp += eigval2[j]*ip*ip;
340 sab += eigval1[i]*tmp;
343 samsb2 = sum1+sum2-2*sab;
349 fprintf(stdout, "The overlap of the covariance matrices:\n");
350 fprintf(stdout, " normalized: %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
351 tmp = 1-sab/sqrt(sum1*sum2);
356 fprintf(stdout, " shape: %.3f\n", 1-sqrt(tmp));
360 static void inprod_matrix(const char *matfile, int natoms,
361 int nvec1, int *eignr1, rvec **eigvec1,
362 int nvec2, int *eignr2, rvec **eigvec2,
363 gmx_bool bSelect, int noutvec, int *outvec)
367 int i, x1, y1, x, y, nlevels;
369 real inp, *t_x, *t_y, max;
377 for (y1 = 0; y1 < nx; y1++)
379 if (outvec[y1] < nvec2)
381 t_y[ny] = eignr2[outvec[y1]]+1;
390 for (y = 0; y < ny; y++)
392 t_y[y] = eignr2[y]+1;
396 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
402 for (x1 = 0; x1 < nx; x1++)
413 t_x[x1] = eignr1[x]+1;
414 fprintf(stderr, " %d", eignr1[x]+1);
415 for (y1 = 0; y1 < ny; y1++)
420 while (outvec[y1] >= nvec2)
430 for (i = 0; i < natoms; i++)
432 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
434 mat[x1][y1] = fabs(inp);
435 if (mat[x1][y1] > max)
441 fprintf(stderr, "\n");
442 rlo.r = 1; rlo.g = 1; rlo.b = 1;
443 rhi.r = 0; rhi.g = 0; rhi.b = 0;
445 out = gmx_ffopen(matfile, "w");
446 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
447 nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
451 static void overlap(const char *outfile, int natoms,
453 int nvec2, int *eignr2, rvec **eigvec2,
454 int noutvec, int *outvec,
455 const output_env_t oenv)
461 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
462 for (i = 0; i < noutvec; i++)
464 fprintf(stderr, "%d ", outvec[i]+1);
466 fprintf(stderr, "\n");
468 out = xvgropen(outfile, "Subspace overlap",
469 "Eigenvectors of trajectory 2", "Overlap", oenv);
470 if (output_env_get_print_xvgr_codes(oenv))
472 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
475 for (x = 0; x < nvec2; x++)
477 for (v = 0; v < noutvec; v++)
481 for (i = 0; i < natoms; i++)
483 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
487 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
493 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
494 const char *projfile, const char *twodplotfile,
495 const char *threedplotfile, const char *filterfile, int skip,
496 const char *extremefile, gmx_bool bExtrAll, real extreme,
497 int nextr, t_atoms *atoms, int natoms, atom_id *index,
498 gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
499 real *sqrtm, rvec *xav,
500 int *eignr, rvec **eigvec,
501 int noutvec, int *outvec, gmx_bool bSplit,
502 const output_env_t oenv)
504 FILE *xvgrout = NULL;
505 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
506 t_trxstatus *out = NULL;
508 int noutvec_extr, imin, imax;
513 real t, inp, **inprod = NULL, min = 0, max = 0;
514 char str[STRLEN], str2[STRLEN], **ylabel, *c;
516 gmx_rmpbc_t gpbc = NULL;
522 noutvec_extr = noutvec;
532 snew(inprod, noutvec+1);
536 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
538 for (i = 0; i < noutvec; i++)
540 fprintf(stderr, "%d ", outvec[i]+1);
542 fprintf(stderr, "\n");
543 out = open_trx(filterfile, "w");
548 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
551 gmx_fatal(FARGS, "the number of atoms in your trajectory (%d) is larger than the number of atoms in your structure file (%d)", nat, atoms->nr);
557 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
560 for (i = 0; i < nat; i++)
570 gmx_rmpbc(gpbc, nat, box, xread);
572 if (nframes >= snew_size)
575 for (i = 0; i < noutvec+1; i++)
577 srenew(inprod[i], snew_size);
580 inprod[noutvec][nframes] = t;
581 /* calculate x: a fitted struture of the selected atoms */
584 reset_x(nfit, ifit, nat, NULL, xread, w_rls);
585 do_fit(nat, w_rls, xref, xread);
587 for (i = 0; i < natoms; i++)
589 copy_rvec(xread[index[i]], x[i]);
592 for (v = 0; v < noutvec; v++)
595 /* calculate (mass-weighted) projection */
597 for (i = 0; i < natoms; i++)
599 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
600 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
601 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
603 inprod[v][nframes] = inp;
607 for (i = 0; i < natoms; i++)
609 for (d = 0; d < DIM; d++)
611 /* misuse xread for output */
612 xread[index[i]][d] = xav[i][d];
613 for (v = 0; v < noutvec; v++)
615 xread[index[i]][d] +=
616 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
620 write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
626 while (read_next_x(oenv, status, &t, xread, box));
636 snew(xread, atoms->nr);
641 gmx_rmpbc_done(gpbc);
647 snew(ylabel, noutvec);
648 for (v = 0; v < noutvec; v++)
650 sprintf(str, "vec %d", eignr[outvec[v]]+1);
651 ylabel[v] = gmx_strdup(str);
653 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
654 write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
655 (const char **)ylabel,
656 nframes, inprod[noutvec], inprod, NULL,
657 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
662 sprintf(str, "projection on eigenvector %d (%s)",
663 eignr[outvec[0]]+1, proj_unit);
664 sprintf(str2, "projection on eigenvector %d (%s)",
665 eignr[outvec[noutvec-1]]+1, proj_unit);
666 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
668 for (i = 0; i < nframes; i++)
670 if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
672 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
674 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
676 gmx_ffclose(xvgrout);
691 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
695 bPDB = fn2ftp(threedplotfile) == efPDB;
697 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
699 b4D = bPDB && (noutvec >= 4);
702 fprintf(stderr, "You have selected four or more eigenvectors:\n"
703 "fourth eigenvector will be plotted "
704 "in bfactor field of pdb file\n");
705 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
706 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
707 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
711 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
712 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
714 init_t_atoms(&atoms, nframes, FALSE);
717 atnm = gmx_strdup("C");
718 resnm = gmx_strdup("PRJ");
722 fact = 10000.0/nframes;
729 for (i = 0; i < nframes; i++)
731 atoms.atomname[i] = &atnm;
732 atoms.atom[i].resind = i;
733 atoms.resinfo[i].name = &resnm;
734 atoms.resinfo[i].nr = ceil(i*fact);
735 atoms.resinfo[i].ic = ' ';
736 x[i][XX] = inprod[0][i];
737 x[i][YY] = inprod[1][i];
738 x[i][ZZ] = inprod[2][i];
744 if ( ( b4D || bSplit ) && bPDB)
746 out = gmx_ffopen(threedplotfile, "w");
747 fprintf(out, "HEADER %s\n", str);
750 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
753 for (i = 0; i < atoms.nr; i++)
755 if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
757 fprintf(out, "TER\n");
760 gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
761 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
764 fprintf(out, "CONECT%5d%5d\n", i, i+1);
768 fprintf(out, "TER\n");
773 write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
775 free_t_atoms(&atoms, FALSE);
780 snew(pmin, noutvec_extr);
781 snew(pmax, noutvec_extr);
784 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
786 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
789 for (v = 0; v < noutvec_extr; v++)
791 for (i = 0; i < nframes; i++)
793 if (inprod[v][i] < inprod[v][imin])
797 if (inprod[v][i] > inprod[v][imax])
802 pmin[v] = inprod[v][imin];
803 pmax[v] = inprod[v][imax];
804 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
806 pmin[v], imin, pmax[v], imax);
814 /* build format string for filename: */
815 strcpy(str, extremefile); /* copy filename */
816 c = strrchr(str, '.'); /* find where extention begins */
817 strcpy(str2, c); /* get extention */
818 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
819 for (v = 0; v < noutvec_extr; v++)
821 /* make filename using format string */
822 if (noutvec_extr == 1)
824 strcpy(str2, extremefile);
828 sprintf(str2, str, eignr[outvec[v]]+1);
830 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
831 nextr, outvec[v]+1, str2);
832 out = open_trx(str2, "w");
833 for (frame = 0; frame < nextr; frame++)
835 if ((extreme == 0) && (nextr <= 3))
837 for (i = 0; i < natoms; i++)
839 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
842 for (i = 0; i < natoms; i++)
844 for (d = 0; d < DIM; d++)
847 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
848 *eigvec[outvec[v]][i][d]/sqrtm[i]);
851 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
858 fprintf(stderr, "\n");
861 static void components(const char *outfile, int natoms,
862 int *eignr, rvec **eigvec,
863 int noutvec, int *outvec,
864 const output_env_t oenv)
868 char str[STRLEN], **ylabel;
870 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
872 snew(ylabel, noutvec);
875 for (i = 0; i < natoms; i++)
879 for (g = 0; g < noutvec; g++)
882 sprintf(str, "vec %d", eignr[v]+1);
883 ylabel[g] = gmx_strdup(str);
885 for (s = 0; s < 4; s++)
887 snew(y[g][s], natoms);
889 for (i = 0; i < natoms; i++)
891 y[g][0][i] = norm(eigvec[v][i]);
892 for (s = 0; s < 3; s++)
894 y[g][s+1][i] = eigvec[v][i][s];
898 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
899 "black: total, red: x, green: y, blue: z",
900 "Atom number", (const char **)ylabel,
901 natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
902 fprintf(stderr, "\n");
905 static void rmsf(const char *outfile, int natoms, real *sqrtm,
906 int *eignr, rvec **eigvec,
907 int noutvec, int *outvec,
908 real *eigval, int neig,
909 const output_env_t oenv)
913 char str[STRLEN], **ylabel;
915 for (i = 0; i < neig; i++)
923 fprintf(stderr, "Writing rmsf to %s\n", outfile);
925 snew(ylabel, noutvec);
928 for (i = 0; i < natoms; i++)
932 for (g = 0; g < noutvec; g++)
935 if (eignr[v] >= neig)
937 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
939 sprintf(str, "vec %d", eignr[v]+1);
940 ylabel[g] = gmx_strdup(str);
942 for (i = 0; i < natoms; i++)
944 y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
947 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
948 "Atom number", (const char **)ylabel,
949 natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
950 fprintf(stderr, "\n");
953 int gmx_anaeig(int argc, char *argv[])
955 static const char *desc[] = {
956 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
957 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
958 "([gmx-nmeig]).[PAR]",
960 "When a trajectory is projected on eigenvectors, all structures are",
961 "fitted to the structure in the eigenvector file, if present, otherwise",
962 "to the structure in the structure file. When no run input file is",
963 "supplied, periodicity will not be taken into account. Most analyses",
964 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
965 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
967 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
968 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
970 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
971 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
973 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
974 "[TT]-first[tt] to [TT]-last[tt].",
975 "The projections of a trajectory on the eigenvectors of its",
976 "covariance matrix are called principal components (pc's).",
977 "It is often useful to check the cosine content of the pc's,",
978 "since the pc's of random diffusion are cosines with the number",
979 "of periods equal to half the pc index.",
980 "The cosine content of the pc's can be calculated with the program",
981 "[gmx-analyze].[PAR]",
983 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
984 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
986 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
987 "three selected eigenvectors.[PAR]",
989 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
990 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
992 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
993 "on the average structure and interpolate [TT]-nframes[tt] frames",
994 "between them, or set your own extremes with [TT]-max[tt]. The",
995 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
996 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
997 "will be written to separate files. Chain identifiers will be added",
998 "when writing a [TT].pdb[tt] file with two or three structures (you",
999 "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1001 " Overlap calculations between covariance analysis:[BR]",
1002 " [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1004 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1005 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1006 "in file [TT]-v[tt].[PAR]",
1008 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1009 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1010 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1011 "have been set explicitly.[PAR]",
1013 "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1014 "a single number for the overlap between the covariance matrices is",
1015 "generated. The formulas are:[BR]",
1016 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1017 "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1018 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1019 "where M1 and M2 are the two covariance matrices and tr is the trace",
1020 "of a matrix. The numbers are proportional to the overlap of the square",
1021 "root of the fluctuations. The normalized overlap is the most useful",
1022 "number, it is 1 for identical matrices and 0 when the sampled",
1023 "subspaces are orthogonal.[PAR]",
1024 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1025 "computed based on the Quasiharmonic approach and based on",
1026 "Schlitter's formula."
1028 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1029 static real max = 0.0, temp = 298.15;
1030 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1032 { "-first", FALSE, etINT, {&first},
1033 "First eigenvector for analysis (-1 is select)" },
1034 { "-last", FALSE, etINT, {&last},
1035 "Last eigenvector for analysis (-1 is till the last)" },
1036 { "-skip", FALSE, etINT, {&skip},
1037 "Only analyse every nr-th frame" },
1038 { "-max", FALSE, etREAL, {&max},
1039 "Maximum for projection of the eigenvector on the average structure, "
1040 "max=0 gives the extremes" },
1041 { "-nframes", FALSE, etINT, {&nextr},
1042 "Number of frames for the extremes output" },
1043 { "-split", FALSE, etBOOL, {&bSplit},
1044 "Split eigenvector projections where time is zero" },
1045 { "-entropy", FALSE, etBOOL, {&bEntropy},
1046 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1047 { "-temp", FALSE, etREAL, {&temp},
1048 "Temperature for entropy calculations" },
1049 { "-nevskip", FALSE, etINT, {&nskip},
1050 "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." }
1052 #define NPA asize(pa)
1058 t_atoms *atoms = NULL;
1059 rvec *xtop, *xref1, *xref2, *xrefp = NULL;
1060 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1061 int nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1062 rvec *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1064 real xid, totmass, *sqrtm, *w_rls, t, lambda;
1067 const char *indexfile;
1070 int nout, *iout, noutvec, *outvec, nfit;
1071 atom_id *index, *ifit;
1072 const char *VecFile, *Vec2File, *topfile;
1073 const char *EigFile, *Eig2File;
1074 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1075 const char *TwoDPlotFile, *ThreeDPlotFile;
1076 const char *FilterFile, *ExtremeFile;
1077 const char *OverlapFile, *InpMatFile;
1078 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1079 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1080 real *eigval1 = NULL, *eigval2 = NULL;
1087 { efTRN, "-v", "eigenvec", ffREAD },
1088 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1089 { efTRX, "-f", NULL, ffOPTRD },
1090 { efTPS, NULL, NULL, ffOPTRD },
1091 { efNDX, NULL, NULL, ffOPTRD },
1092 { efXVG, "-eig", "eigenval", ffOPTRD },
1093 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1094 { efXVG, "-comp", "eigcomp", ffOPTWR },
1095 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1096 { efXVG, "-proj", "proj", ffOPTWR },
1097 { efXVG, "-2d", "2dproj", ffOPTWR },
1098 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1099 { efTRX, "-filt", "filtered", ffOPTWR },
1100 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1101 { efXVG, "-over", "overlap", ffOPTWR },
1102 { efXPM, "-inpr", "inprod", ffOPTWR }
1104 #define NFILE asize(fnm)
1106 if (!parse_common_args(&argc, argv,
1107 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
1108 NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1113 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1115 VecFile = opt2fn("-v", NFILE, fnm);
1116 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1117 topfile = ftp2fn(efTPS, NFILE, fnm);
1118 EigFile = opt2fn_null("-eig", NFILE, fnm);
1119 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1120 CompFile = opt2fn_null("-comp", NFILE, fnm);
1121 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1122 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1123 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1124 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1125 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1126 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1127 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1128 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1130 bTop = fn2bTPX(topfile);
1131 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1132 || FilterFile || ExtremeFile;
1134 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1135 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1136 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1137 bVec2 = Vec2File || OverlapFile || InpMatFile;
1138 bM = RmsfFile || bProj;
1139 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1140 || TwoDPlotFile || ThreeDPlotFile;
1141 bIndex = bM || bProj;
1142 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1143 FilterFile || (bIndex && indexfile);
1144 bCompare = Vec2File || Eig2File;
1145 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1147 read_eigenvectors(VecFile, &natoms, &bFit1,
1148 &xref1, &bDMR1, &xav1, &bDMA1,
1149 &nvec1, &eignr1, &eigvec1, &eigval1);
1152 /* Overwrite eigenvalues from separate files if the user provides them */
1153 if (EigFile != NULL)
1155 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1156 if (neig_tmp != neig1)
1158 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1161 srenew(eigval1, neig1);
1162 for (j = 0; j < neig1; j++)
1164 real tmp = eigval1[j];
1165 eigval1[j] = xvgdata[1][j];
1166 if (debug && (eigval1[j] != tmp))
1168 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1169 j, tmp, eigval1[j]);
1172 for (j = 0; j < i; j++)
1177 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1184 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1186 calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1187 calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1194 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1196 read_eigenvectors(Vec2File, &neig2, &bFit2,
1197 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1202 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1206 if (Eig2File != NULL)
1208 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1209 srenew(eigval2, neig2);
1210 for (j = 0; j < neig2; j++)
1212 eigval2[j] = xvgdata[1][j];
1214 for (j = 0; j < i; j++)
1219 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1223 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1227 if ((xref1 == NULL) && (bM || bTraj))
1243 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1244 title, &top, &ePBC, &xtop, NULL, topbox, bM);
1246 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1247 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1248 /* Fitting is only required for the projection */
1253 printf("\nNote: the structure in %s should be the same\n"
1254 " as the one used for the fit in g_covar\n", topfile);
1256 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1257 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1259 snew(w_rls, atoms->nr);
1260 for (i = 0; (i < nfit); i++)
1264 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1268 w_rls[ifit[i]] = 1.0;
1272 snew(xrefp, atoms->nr);
1275 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1278 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);
1280 for (i = 0; (i < nfit); i++)
1282 copy_rvec(xref1[i], xrefp[ifit[i]]);
1287 /* The top coordinates are the fitting reference */
1288 for (i = 0; (i < nfit); i++)
1290 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1292 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1295 gmx_rmpbc_done(gpbc);
1300 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1301 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1304 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1309 snew(sqrtm, natoms);
1312 proj_unit = "u\\S1/2\\Nnm";
1313 for (i = 0; (i < natoms); i++)
1315 sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1321 for (i = 0; (i < natoms); i++)
1331 for (i = 0; (i < natoms); i++)
1333 for (d = 0; (d < DIM); d++)
1335 t += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1336 totmass += sqr(sqrtm[i]);
1339 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1340 " %.3f (nm)\n\n", sqrt(t/totmass));
1351 /* make an index from first to last */
1352 nout = last-first+1;
1354 for (i = 0; i < nout; i++)
1356 iout[i] = first-1+i;
1359 else if (ThreeDPlotFile)
1361 /* make an index of first+(0,1,2) and last */
1362 nout = bPDB3D ? 4 : 3;
1363 nout = min(last-first+1, nout);
1371 iout[nout-1] = last-1;
1375 /* make an index of first and last */
1384 printf("Select eigenvectors for output, end your selection with 0\n");
1391 srenew(iout, nout+1);
1392 if (1 != scanf("%d", &iout[nout]))
1394 gmx_fatal(FARGS, "Error reading user input");
1398 while (iout[nout] >= 0);
1402 /* make an index of the eigenvectors which are present */
1405 for (i = 0; i < nout; i++)
1408 while ((j < nvec1) && (eignr1[j] != iout[i]))
1412 if ((j < nvec1) && (eignr1[j] == iout[i]))
1414 outvec[noutvec] = j;
1418 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1421 fprintf(stderr, ":");
1422 for (j = 0; j < noutvec; j++)
1424 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1427 fprintf(stderr, "\n");
1431 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1436 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1442 project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1443 bTop ? &top : NULL, ePBC, topbox,
1444 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1445 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1446 bFit1, xrefp, nfit, ifit, w_rls,
1447 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1453 overlap(OverlapFile, natoms,
1454 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1459 inprod_matrix(InpMatFile, natoms,
1460 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1461 bFirstLastSet, noutvec, outvec);
1466 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1470 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1471 !bCompare && !bEntropy)
1473 fprintf(stderr, "\nIf you want some output,"
1474 " set one (or two or ...) of the output file options\n");
1478 view_all(oenv, NFILE, fnm);