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,2015,2016,2017,2018, 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.
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/commandline/viewit.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/fileio/matio.h"
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/fileio/tpxio.h"
53 #include "gromacs/fileio/trxio.h"
54 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/gmxana/eigio.h"
56 #include "gromacs/gmxana/gmx_ana.h"
57 #include "gromacs/math/do_fit.h"
58 #include "gromacs/math/functions.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/pbcutil/rmpbc.h"
61 #include "gromacs/topology/index.h"
62 #include "gromacs/topology/topology.h"
63 #include "gromacs/utility/arraysize.h"
64 #include "gromacs/utility/cstringutil.h"
65 #include "gromacs/utility/fatalerror.h"
66 #include "gromacs/utility/futil.h"
67 #include "gromacs/utility/gmxassert.h"
68 #include "gromacs/utility/smalloc.h"
70 #include "thermochemistry.h"
72 static const char *proj_unit;
74 static real tick_spacing(real range, int minticks)
83 sp = 0.2*std::exp(std::log(10.0)*std::ceil(std::log(range)/std::log(10.0)));
84 while (range/sp < minticks-1)
92 static void write_xvgr_graphs(const char *file, int ngraphs, int nsetspergraph,
93 const char *title, const char *subtitle,
94 const std::string &xlabel, const char **ylabel,
95 int n, real *x, real **y, real ***sy,
96 real scale_x, gmx_bool bZero, gmx_bool bSplit,
97 const gmx_output_env_t *oenv)
101 real ymin, ymax, xsp, ysp;
103 out = gmx_ffopen(file, "w");
104 if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
106 fprintf(out, "@ autoscale onread none\n");
108 for (g = 0; g < ngraphs; g++)
114 for (i = 0; i < n; i++)
131 for (s = 0; s < nsetspergraph; s++)
133 for (i = 0; i < n; i++)
135 if (sy[g][s][i] < ymin)
139 if (sy[g][s][i] > ymax)
152 ymin = ymin-0.1*(ymax-ymin);
154 ymax = ymax+0.1*(ymax-ymin);
155 xsp = tick_spacing((x[n-1]-x[0])*scale_x, 4);
156 ysp = tick_spacing(ymax-ymin, 3);
157 if (output_env_get_print_xvgr_codes(oenv))
159 fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
162 fprintf(out, "@ title \"%s\"\n", title);
165 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
170 fprintf(out, "@ xaxis label \"%s\"\n", xlabel.c_str());
174 fprintf(out, "@ xaxis ticklabel off\n");
178 fprintf(out, "@ world xmin %g\n", x[0]*scale_x);
179 fprintf(out, "@ world xmax %g\n", x[n-1]*scale_x);
180 fprintf(out, "@ world ymin %g\n", ymin);
181 fprintf(out, "@ world ymax %g\n", ymax);
183 fprintf(out, "@ view xmin 0.15\n");
184 fprintf(out, "@ view xmax 0.85\n");
185 fprintf(out, "@ view ymin %g\n", 0.15+(ngraphs-1-g)*0.7/ngraphs);
186 fprintf(out, "@ view ymax %g\n", 0.15+(ngraphs-g)*0.7/ngraphs);
187 fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
188 fprintf(out, "@ xaxis tick major %g\n", xsp);
189 fprintf(out, "@ xaxis tick minor %g\n", xsp/2);
190 fprintf(out, "@ xaxis ticklabel start type spec\n");
191 fprintf(out, "@ xaxis ticklabel start %g\n", std::ceil(ymin/xsp)*xsp);
192 fprintf(out, "@ yaxis tick major %g\n", ysp);
193 fprintf(out, "@ yaxis tick minor %g\n", ysp/2);
194 fprintf(out, "@ yaxis ticklabel start type spec\n");
195 fprintf(out, "@ yaxis ticklabel start %g\n", std::ceil(ymin/ysp)*ysp);
196 if ((ymin < 0) && (ymax > 0))
198 fprintf(out, "@ zeroxaxis bar on\n");
199 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
202 for (s = 0; s < nsetspergraph; s++)
204 for (i = 0; i < n; i++)
206 if (bSplit && i > 0 && std::abs(x[i]) < 1e-5)
208 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
210 fprintf(out, "%10.4f %10.5f\n",
211 x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
213 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
220 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
221 real *eigval1, int neig1, real *eigval2, int neig2)
225 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
227 n = std::min(n1, n2);
229 n = std::min(n, std::min(neig1, neig2));
230 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
233 for (i = 0; i < n; i++)
240 eigval1[i] = std::sqrt(eigval1[i]);
243 for (i = n; i < neig1; i++)
245 trace1 += eigval1[i];
248 for (i = 0; i < n; i++)
255 eigval2[i] = std::sqrt(eigval2[i]);
260 // The static analyzer appears to be confused by the fact that the loop below
261 // starts from n instead of 0. However, given all the complex code it's
262 // better to be safe than sorry, so we check it with an assert.
263 // If we are in this comparison routine in the first place, neig2 should not be 0,
264 // so eigval2 should always be a valid pointer.
265 GMX_RELEASE_ASSERT(eigval2 != nullptr, "NULL pointer provided for eigval2");
267 for (i = n; i < neig2; i++)
269 trace2 += eigval2[i];
272 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
273 if (neig1 != n || neig2 != n)
275 fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
276 static_cast<int>(100*sum1/trace1+0.5), static_cast<int>(100*sum2/trace2+0.5));
278 fprintf(stdout, "Square root of the traces: %g and %g\n",
279 std::sqrt(sum1), std::sqrt(sum2));
282 for (i = 0; i < n; i++)
285 for (j = 0; j < n; j++)
288 for (k = 0; k < natoms; k++)
290 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
292 tmp += eigval2[j]*ip*ip;
294 sab += eigval1[i]*tmp;
297 samsb2 = sum1+sum2-2*sab;
303 fprintf(stdout, "The overlap of the covariance matrices:\n");
304 fprintf(stdout, " normalized: %.3f\n", 1-std::sqrt(samsb2/(sum1+sum2)));
305 tmp = 1-sab/std::sqrt(sum1*sum2);
310 fprintf(stdout, " shape: %.3f\n", 1-std::sqrt(tmp));
314 static void inprod_matrix(const char *matfile, int natoms,
315 int nvec1, int *eignr1, rvec **eigvec1,
316 int nvec2, const int *eignr2, rvec **eigvec2,
317 gmx_bool bSelect, int noutvec, const int *outvec)
321 int i, x1, y1, x, y, nlevels;
323 real inp, *t_x, *t_y, maxval;
331 for (y1 = 0; y1 < nx; y1++)
333 if (outvec[y1] < nvec2)
335 t_y[ny] = eignr2[outvec[y1]]+1;
344 for (y = 0; y < ny; y++)
346 t_y[y] = eignr2[y]+1;
350 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
356 for (x1 = 0; x1 < nx; x1++)
367 t_x[x1] = eignr1[x]+1;
368 fprintf(stderr, " %d", eignr1[x]+1);
369 for (y1 = 0; y1 < ny; y1++)
374 while (outvec[y1] >= nvec2)
384 for (i = 0; i < natoms; i++)
386 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
388 mat[x1][y1] = std::abs(inp);
389 if (mat[x1][y1] > maxval)
391 maxval = mat[x1][y1];
395 fprintf(stderr, "\n");
396 rlo.r = 1; rlo.g = 1; rlo.b = 1;
397 rhi.r = 0; rhi.g = 0; rhi.b = 0;
399 out = gmx_ffopen(matfile, "w");
400 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
401 nx, ny, t_x, t_y, mat, 0.0, maxval, rlo, rhi, &nlevels);
405 static void overlap(const char *outfile, int natoms,
407 int nvec2, int *eignr2, rvec **eigvec2,
408 int noutvec, int *outvec,
409 const gmx_output_env_t *oenv)
415 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
416 for (i = 0; i < noutvec; i++)
418 fprintf(stderr, "%d ", outvec[i]+1);
420 fprintf(stderr, "\n");
422 out = xvgropen(outfile, "Subspace overlap",
423 "Eigenvectors of trajectory 2", "Overlap", oenv);
424 if (output_env_get_print_xvgr_codes(oenv))
426 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
429 for (x = 0; x < nvec2; x++)
431 for (v = 0; v < noutvec; v++)
435 for (i = 0; i < natoms; i++)
437 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
439 overlap += gmx::square(inp);
441 fprintf(out, "%5d %5.3f\n", eignr2[x]+1, overlap/noutvec);
447 static void project(const char *trajfile, const t_topology *top, int ePBC, matrix topbox,
448 const char *projfile, const char *twodplotfile,
449 const char *threedplotfile, const char *filterfile, int skip,
450 const char *extremefile, gmx_bool bExtrAll, real extreme,
451 int nextr, const t_atoms *atoms, int natoms, int *index,
452 gmx_bool bFit, rvec *xref, int nfit, int *ifit, real *w_rls,
453 const real *sqrtm, rvec *xav,
454 int *eignr, rvec **eigvec,
455 int noutvec, int *outvec, gmx_bool bSplit,
456 const gmx_output_env_t *oenv)
458 FILE *xvgrout = nullptr;
459 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
460 t_trxstatus *out = nullptr;
462 int noutvec_extr, imin, imax;
467 real t, inp, **inprod = nullptr;
468 char str[STRLEN], str2[STRLEN], *c;
471 gmx_rmpbc_t gpbc = nullptr;
477 noutvec_extr = noutvec;
487 snew(inprod, noutvec+1);
491 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
493 for (i = 0; i < noutvec; i++)
495 fprintf(stderr, "%d ", outvec[i]+1);
497 fprintf(stderr, "\n");
498 out = open_trx(filterfile, "w");
503 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
506 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);
512 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
515 for (i = 0; i < nat; i++)
525 gmx_rmpbc(gpbc, nat, box, xread);
527 if (nframes >= snew_size)
530 for (i = 0; i < noutvec+1; i++)
532 srenew(inprod[i], snew_size);
535 inprod[noutvec][nframes] = t;
536 /* calculate x: a fitted struture of the selected atoms */
539 reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
540 do_fit(nat, w_rls, xref, xread);
542 for (i = 0; i < natoms; i++)
544 copy_rvec(xread[index[i]], x[i]);
547 for (v = 0; v < noutvec; v++)
550 /* calculate (mass-weighted) projection */
552 for (i = 0; i < natoms; i++)
554 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
555 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
556 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
558 inprod[v][nframes] = inp;
562 for (i = 0; i < natoms; i++)
564 for (d = 0; d < DIM; d++)
566 /* misuse xread for output */
567 xread[index[i]][d] = xav[i][d];
568 for (v = 0; v < noutvec; v++)
570 xread[index[i]][d] +=
571 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
575 write_trx(out, natoms, index, atoms, 0, t, box, xread, nullptr, nullptr);
581 while (read_next_x(oenv, status, &t, xread, box));
591 snew(xread, atoms->nr);
596 gmx_rmpbc_done(gpbc);
602 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL if projfile is non-NULL");
603 snew(ylabel, noutvec);
604 for (v = 0; v < noutvec; v++)
606 sprintf(str, "vec %d", eignr[outvec[v]]+1);
607 ylabel[v] = gmx_strdup(str);
609 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
610 write_xvgr_graphs(projfile, noutvec, 1, str, nullptr, output_env_get_xvgr_tlabel(oenv),
612 nframes, inprod[noutvec], inprod, nullptr,
613 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
618 sprintf(str, "projection on eigenvector %d (%s)",
619 eignr[outvec[0]]+1, proj_unit);
620 sprintf(str2, "projection on eigenvector %d (%s)",
621 eignr[outvec[noutvec-1]]+1, proj_unit);
622 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
624 for (i = 0; i < nframes; i++)
626 if (bSplit && i > 0 && std::abs(inprod[noutvec][i]) < 1e-5)
628 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
630 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
647 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
651 bPDB = fn2ftp(threedplotfile) == efPDB;
653 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
655 b4D = bPDB && (noutvec >= 4);
658 fprintf(stderr, "You have selected four or more eigenvectors:\n"
659 "fourth eigenvector will be plotted "
660 "in bfactor field of pdb file\n");
661 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
662 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
663 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
667 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
668 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
670 init_t_atoms(&atoms, nframes, FALSE);
673 atnm = gmx_strdup("C");
674 resnm = gmx_strdup("PRJ");
678 fact = 10000.0/nframes;
685 for (i = 0; i < nframes; i++)
687 atoms.atomname[i] = &atnm;
688 atoms.atom[i].resind = i;
689 atoms.resinfo[i].name = &resnm;
690 atoms.resinfo[i].nr = static_cast<int>(std::ceil(i*fact));
691 atoms.resinfo[i].ic = ' ';
692 x[i][XX] = inprod[0][i];
693 x[i][YY] = inprod[1][i];
694 x[i][ZZ] = inprod[2][i];
700 if ( ( b4D || bSplit ) && bPDB)
702 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL with 4D or split PDB output options");
704 out = gmx_ffopen(threedplotfile, "w");
705 fprintf(out, "HEADER %s\n", str);
708 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
711 for (i = 0; i < atoms.nr; i++)
713 if (j > 0 && bSplit && std::abs(inprod[noutvec][i]) < 1e-5)
715 fprintf(out, "TER\n");
718 gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
719 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
722 fprintf(out, "CONECT%5d%5d\n", i, i+1);
726 fprintf(out, "TER\n");
731 write_sto_conf(threedplotfile, str, &atoms, x, nullptr, ePBC, box);
738 snew(pmin, noutvec_extr);
739 snew(pmax, noutvec_extr);
742 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL");
743 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
745 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
748 for (v = 0; v < noutvec_extr; v++)
750 for (i = 0; i < nframes; i++)
752 if (inprod[v][i] < inprod[v][imin])
756 if (inprod[v][i] > inprod[v][imax])
761 pmin[v] = inprod[v][imin];
762 pmax[v] = inprod[v][imax];
763 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
765 pmin[v], imin, pmax[v], imax);
773 /* build format string for filename: */
774 std::strcpy(str, extremefile); /* copy filename */
775 c = std::strrchr(str, '.'); /* find where extention begins */
776 std::strcpy(str2, c); /* get extention */
777 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
778 for (v = 0; v < noutvec_extr; v++)
780 /* make filename using format string */
781 if (noutvec_extr == 1)
783 std::strcpy(str2, extremefile);
787 sprintf(str2, str, eignr[outvec[v]]+1);
789 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
790 nextr, outvec[v]+1, str2);
791 out = open_trx(str2, "w");
792 for (frame = 0; frame < nextr; frame++)
794 if ((extreme == 0) && (nextr <= 3))
796 for (i = 0; i < natoms; i++)
798 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
801 for (i = 0; i < natoms; i++)
803 for (d = 0; d < DIM; d++)
806 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
807 *eigvec[outvec[v]][i][d]/sqrtm[i]);
810 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, nullptr, nullptr);
817 fprintf(stderr, "\n");
820 static void components(const char *outfile, int natoms,
821 int *eignr, rvec **eigvec,
822 int noutvec, const int *outvec,
823 const gmx_output_env_t *oenv)
830 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
832 snew(ylabel, noutvec);
835 for (i = 0; i < natoms; i++)
839 for (g = 0; g < noutvec; g++)
842 sprintf(str, "vec %d", eignr[v]+1);
843 ylabel[g] = gmx_strdup(str);
845 for (s = 0; s < 4; s++)
847 snew(y[g][s], natoms);
849 for (i = 0; i < natoms; i++)
851 y[g][0][i] = norm(eigvec[v][i]);
852 for (s = 0; s < 3; s++)
854 y[g][s+1][i] = eigvec[v][i][s];
858 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
859 "black: total, red: x, green: y, blue: z",
860 "Atom number", ylabel,
861 natoms, x, nullptr, y, 1, FALSE, FALSE, oenv);
862 fprintf(stderr, "\n");
865 static void rmsf(const char *outfile, int natoms, const real *sqrtm,
866 int *eignr, rvec **eigvec,
867 int noutvec, const int *outvec,
868 real *eigval, int neig,
869 const gmx_output_env_t *oenv)
876 for (i = 0; i < neig; i++)
884 fprintf(stderr, "Writing rmsf to %s\n", outfile);
886 snew(ylabel, noutvec);
889 for (i = 0; i < natoms; i++)
893 for (g = 0; g < noutvec; g++)
896 if (eignr[v] >= neig)
898 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
900 sprintf(str, "vec %d", eignr[v]+1);
901 ylabel[g] = gmx_strdup(str);
903 for (i = 0; i < natoms; i++)
905 y[g][i] = std::sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
908 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", nullptr,
909 "Atom number", ylabel,
910 natoms, x, y, nullptr, 1, TRUE, FALSE, oenv);
911 fprintf(stderr, "\n");
914 int gmx_anaeig(int argc, char *argv[])
916 static const char *desc[] = {
917 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
918 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
919 "([gmx-nmeig]).[PAR]",
921 "When a trajectory is projected on eigenvectors, all structures are",
922 "fitted to the structure in the eigenvector file, if present, otherwise",
923 "to the structure in the structure file. When no run input file is",
924 "supplied, periodicity will not be taken into account. Most analyses",
925 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
926 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
928 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
929 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
931 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
932 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
934 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
935 "[TT]-first[tt] to [TT]-last[tt].",
936 "The projections of a trajectory on the eigenvectors of its",
937 "covariance matrix are called principal components (pc's).",
938 "It is often useful to check the cosine content of the pc's,",
939 "since the pc's of random diffusion are cosines with the number",
940 "of periods equal to half the pc index.",
941 "The cosine content of the pc's can be calculated with the program",
942 "[gmx-analyze].[PAR]",
944 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
945 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
947 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
948 "three selected eigenvectors.[PAR]",
950 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
951 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
953 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
954 "on the average structure and interpolate [TT]-nframes[tt] frames",
955 "between them, or set your own extremes with [TT]-max[tt]. The",
956 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
957 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
958 "will be written to separate files. Chain identifiers will be added",
959 "when writing a [REF].pdb[ref] file with two or three structures (you",
960 "can use [TT]rasmol -nmrpdb[tt] to view such a [REF].pdb[ref] file).[PAR]",
962 "Overlap calculations between covariance analysis",
963 "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^",
965 "[BB]Note:[bb] the analysis should use the same fitting structure",
967 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
968 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
969 "in file [TT]-v[tt].[PAR]",
971 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
972 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
973 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
974 "have been set explicitly.[PAR]",
976 "When [TT]-v[tt] and [TT]-v2[tt] are given, a single number for the",
977 "overlap between the covariance matrices is generated. Note that the",
978 "eigenvalues are by default read from the timestamp field in the",
979 "eigenvector input files, but when [TT]-eig[tt], or [TT]-eig2[tt] are",
980 "given, the corresponding eigenvalues are used instead. The formulas are::",
982 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))",
983 " normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))",
984 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))",
986 "where M1 and M2 are the two covariance matrices and tr is the trace",
987 "of a matrix. The numbers are proportional to the overlap of the square",
988 "root of the fluctuations. The normalized overlap is the most useful",
989 "number, it is 1 for identical matrices and 0 when the sampled",
990 "subspaces are orthogonal.[PAR]",
991 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
992 "computed based on the Quasiharmonic approach and based on",
993 "Schlitter's formula."
995 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
996 static real max = 0.0, temp = 298.15;
997 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
999 { "-first", FALSE, etINT, {&first},
1000 "First eigenvector for analysis (-1 is select)" },
1001 { "-last", FALSE, etINT, {&last},
1002 "Last eigenvector for analysis (-1 is till the last)" },
1003 { "-skip", FALSE, etINT, {&skip},
1004 "Only analyse every nr-th frame" },
1005 { "-max", FALSE, etREAL, {&max},
1006 "Maximum for projection of the eigenvector on the average structure, "
1007 "max=0 gives the extremes" },
1008 { "-nframes", FALSE, etINT, {&nextr},
1009 "Number of frames for the extremes output" },
1010 { "-split", FALSE, etBOOL, {&bSplit},
1011 "Split eigenvector projections where time is zero" },
1012 { "-entropy", FALSE, etBOOL, {&bEntropy},
1013 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1014 { "-temp", FALSE, etREAL, {&temp},
1015 "Temperature for entropy calculations" },
1016 { "-nevskip", FALSE, etINT, {&nskip},
1017 "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." }
1019 #define NPA asize(pa)
1023 const t_atoms *atoms = nullptr;
1024 rvec *xtop, *xref1, *xref2, *xrefp = nullptr;
1025 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1026 int nvec1, nvec2, *eignr1 = nullptr, *eignr2 = nullptr;
1027 rvec *xav1, *xav2, **eigvec1 = nullptr, **eigvec2 = nullptr;
1029 real totmass, *sqrtm, *w_rls, t;
1032 const char *indexfile;
1034 int nout, *iout, noutvec, *outvec, nfit;
1035 int *index = nullptr, *ifit = nullptr;
1036 const char *VecFile, *Vec2File, *topfile;
1037 const char *EigFile, *Eig2File;
1038 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1039 const char *TwoDPlotFile, *ThreeDPlotFile;
1040 const char *FilterFile, *ExtremeFile;
1041 const char *OverlapFile, *InpMatFile;
1042 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1043 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1044 real *eigval1 = nullptr, *eigval2 = nullptr;
1047 gmx_output_env_t *oenv;
1051 { efTRN, "-v", "eigenvec", ffREAD },
1052 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1053 { efTRX, "-f", nullptr, ffOPTRD },
1054 { efTPS, nullptr, nullptr, ffOPTRD },
1055 { efNDX, nullptr, nullptr, ffOPTRD },
1056 { efXVG, "-eig", "eigenval", ffOPTRD },
1057 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1058 { efXVG, "-comp", "eigcomp", ffOPTWR },
1059 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1060 { efXVG, "-proj", "proj", ffOPTWR },
1061 { efXVG, "-2d", "2dproj", ffOPTWR },
1062 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1063 { efTRX, "-filt", "filtered", ffOPTWR },
1064 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1065 { efXVG, "-over", "overlap", ffOPTWR },
1066 { efXPM, "-inpr", "inprod", ffOPTWR }
1068 #define NFILE asize(fnm)
1070 if (!parse_common_args(&argc, argv,
1071 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
1072 NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1077 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1079 VecFile = opt2fn("-v", NFILE, fnm);
1080 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1081 topfile = ftp2fn(efTPS, NFILE, fnm);
1082 EigFile = opt2fn_null("-eig", NFILE, fnm);
1083 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1084 CompFile = opt2fn_null("-comp", NFILE, fnm);
1085 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1086 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1087 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1088 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1089 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1090 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1091 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1092 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1094 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1095 || FilterFile || ExtremeFile;
1097 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1098 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1099 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1100 bVec2 = Vec2File || OverlapFile || InpMatFile;
1101 bM = RmsfFile || bProj;
1102 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1103 || TwoDPlotFile || ThreeDPlotFile;
1104 bIndex = bM || bProj;
1105 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1106 FilterFile || (bIndex && indexfile);
1107 bCompare = Vec2File || Eig2File;
1108 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1110 read_eigenvectors(VecFile, &natoms, &bFit1,
1111 &xref1, &bDMR1, &xav1, &bDMA1,
1112 &nvec1, &eignr1, &eigvec1, &eigval1);
1115 /* Overwrite eigenvalues from separate files if the user provides them */
1116 if (EigFile != nullptr)
1118 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1119 if (neig_tmp != neig1)
1121 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1124 srenew(eigval1, neig1);
1125 for (j = 0; j < neig1; j++)
1127 real tmp = eigval1[j];
1128 eigval1[j] = xvgdata[1][j];
1129 if (debug && (eigval1[j] != tmp))
1131 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1132 j, tmp, eigval1[j]);
1135 for (j = 0; j < i; j++)
1140 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1147 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1149 printf("The Entropy due to the Schlitter formula is %g J/mol K\n",
1150 calcSchlitterEntropy(gmx::arrayRefFromArray(eigval1, neig1),
1158 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1160 read_eigenvectors(Vec2File, &neig2, &bFit2,
1161 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1166 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1175 if (Eig2File != nullptr)
1177 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1178 srenew(eigval2, neig2);
1179 for (j = 0; j < neig2; j++)
1181 eigval2[j] = xvgdata[1][j];
1183 for (j = 0; j < i; j++)
1188 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1192 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1196 if ((xref1 == nullptr) && (bM || bTraj))
1212 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1213 &top, &ePBC, &xtop, nullptr, topbox, bM);
1215 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1216 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1217 /* Fitting is only required for the projection */
1220 if (xref1 == nullptr)
1222 printf("\nNote: the structure in %s should be the same\n"
1223 " as the one used for the fit in g_covar\n", topfile);
1225 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1226 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1228 snew(w_rls, atoms->nr);
1229 for (i = 0; (i < nfit); i++)
1233 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1237 w_rls[ifit[i]] = 1.0;
1241 snew(xrefp, atoms->nr);
1242 if (xref1 != nullptr)
1244 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1247 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);
1249 for (i = 0; (i < nfit); i++)
1251 copy_rvec(xref1[i], xrefp[ifit[i]]);
1256 /* The top coordinates are the fitting reference */
1257 for (i = 0; (i < nfit); i++)
1259 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1261 reset_x(nfit, ifit, atoms->nr, nullptr, xrefp, w_rls);
1264 gmx_rmpbc_done(gpbc);
1269 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1270 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1273 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1278 snew(sqrtm, natoms);
1281 proj_unit = "u\\S1/2\\Nnm";
1282 for (i = 0; (i < natoms); i++)
1284 sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
1290 for (i = 0; (i < natoms); i++)
1300 for (i = 0; (i < natoms); i++)
1302 for (d = 0; (d < DIM); d++)
1304 t += gmx::square((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1305 totmass += gmx::square(sqrtm[i]);
1308 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1309 " %.3f (nm)\n\n", std::sqrt(t/totmass));
1320 /* make an index from first to last */
1321 nout = last-first+1;
1323 for (i = 0; i < nout; i++)
1325 iout[i] = first-1+i;
1328 else if (ThreeDPlotFile)
1330 /* make an index of first+(0,1,2) and last */
1331 nout = bPDB3D ? 4 : 3;
1332 nout = std::min(last-first+1, nout);
1340 iout[nout-1] = last-1;
1344 /* make an index of first and last */
1353 printf("Select eigenvectors for output, end your selection with 0\n");
1360 srenew(iout, nout+1);
1361 if (1 != scanf("%d", &iout[nout]))
1363 gmx_fatal(FARGS, "Error reading user input");
1367 while (iout[nout] >= 0);
1371 /* make an index of the eigenvectors which are present */
1374 for (i = 0; i < nout; i++)
1377 while ((j < nvec1) && (eignr1[j] != iout[i]))
1381 if ((j < nvec1) && (eignr1[j] == iout[i]))
1383 outvec[noutvec] = j;
1387 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1390 fprintf(stderr, ":");
1391 for (j = 0; j < noutvec; j++)
1393 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1396 fprintf(stderr, "\n");
1400 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1405 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1411 project(bTraj ? opt2fn("-f", NFILE, fnm) : nullptr,
1412 bTop ? &top : nullptr, ePBC, topbox,
1413 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1414 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1415 bFit1, xrefp, nfit, ifit, w_rls,
1416 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1422 overlap(OverlapFile, natoms,
1423 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1428 inprod_matrix(InpMatFile, natoms,
1429 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1430 bFirstLastSet, noutvec, outvec);
1435 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1439 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1440 !bCompare && !bEntropy)
1442 fprintf(stderr, "\nIf you want some output,"
1443 " set one (or two or ...) of the output file options\n");
1447 view_all(oenv, NFILE, fnm);