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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/commandline/viewit.h"
50 #include "gromacs/fileio/confio.h"
51 #include "gromacs/fileio/matio.h"
52 #include "gromacs/fileio/pdbio.h"
53 #include "gromacs/fileio/tpxio.h"
54 #include "gromacs/fileio/trxio.h"
55 #include "gromacs/fileio/xvgr.h"
56 #include "gromacs/gmxana/eigio.h"
57 #include "gromacs/gmxana/gmx_ana.h"
58 #include "gromacs/math/do_fit.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/pbcutil/rmpbc.h"
62 #include "gromacs/topology/index.h"
63 #include "gromacs/topology/topology.h"
64 #include "gromacs/utility/arraysize.h"
65 #include "gromacs/utility/cstringutil.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/futil.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/smalloc.h"
71 #include "thermochemistry.h"
73 static const char* proj_unit;
75 static real tick_spacing(real range, int minticks)
84 sp = 0.2 * std::exp(std::log(10.0) * std::ceil(std::log(range) / std::log(10.0)));
85 while (range / sp < minticks - 1)
93 static void write_xvgr_graphs(const char* file,
98 const std::string& xlabel,
107 const gmx_output_env_t* oenv)
111 real ymin, ymax, xsp, ysp;
113 out = gmx_ffopen(file, "w");
114 if (output_env_get_xvg_format(oenv) == XvgFormat::Xmgrace)
116 fprintf(out, "@ autoscale onread none\n");
118 for (g = 0; g < ngraphs; g++)
124 for (i = 0; i < n; i++)
141 for (s = 0; s < nsetspergraph; s++)
143 for (i = 0; i < n; i++)
145 if (sy[g][s][i] < ymin)
149 if (sy[g][s][i] > ymax)
162 ymin = ymin - 0.1 * (ymax - ymin);
164 ymax = ymax + 0.1 * (ymax - ymin);
165 xsp = tick_spacing((x[n - 1] - x[0]) * scale_x, 4);
166 ysp = tick_spacing(ymax - ymin, 3);
167 if (output_env_get_print_xvgr_codes(oenv))
169 fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
172 fprintf(out, "@ title \"%s\"\n", title);
175 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
178 if (g == ngraphs - 1)
180 fprintf(out, "@ xaxis label \"%s\"\n", xlabel.c_str());
184 fprintf(out, "@ xaxis ticklabel off\n");
188 fprintf(out, "@ world xmin %g\n", x[0] * scale_x);
189 fprintf(out, "@ world xmax %g\n", x[n - 1] * scale_x);
190 fprintf(out, "@ world ymin %g\n", ymin);
191 fprintf(out, "@ world ymax %g\n", ymax);
193 fprintf(out, "@ view xmin 0.15\n");
194 fprintf(out, "@ view xmax 0.85\n");
195 fprintf(out, "@ view ymin %g\n", 0.15 + (ngraphs - 1 - g) * 0.7 / ngraphs);
196 fprintf(out, "@ view ymax %g\n", 0.15 + (ngraphs - g) * 0.7 / ngraphs);
197 fprintf(out, "@ yaxis label \"%s\"\n", ylabel[g]);
198 fprintf(out, "@ xaxis tick major %g\n", xsp);
199 fprintf(out, "@ xaxis tick minor %g\n", xsp / 2);
200 fprintf(out, "@ xaxis ticklabel start type spec\n");
201 fprintf(out, "@ xaxis ticklabel start %g\n", std::ceil(ymin / xsp) * xsp);
202 fprintf(out, "@ yaxis tick major %g\n", ysp);
203 fprintf(out, "@ yaxis tick minor %g\n", ysp / 2);
204 fprintf(out, "@ yaxis ticklabel start type spec\n");
205 fprintf(out, "@ yaxis ticklabel start %g\n", std::ceil(ymin / ysp) * ysp);
206 if ((ymin < 0) && (ymax > 0))
208 fprintf(out, "@ zeroxaxis bar on\n");
209 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
212 for (s = 0; s < nsetspergraph; s++)
214 for (i = 0; i < n; i++)
216 if (bSplit && i > 0 && std::abs(x[i]) < 1e-5)
218 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
220 fprintf(out, "%10.4f %10.5f\n", x[i] * scale_x, y ? y[g][i] : sy[g][s][i]);
222 fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
229 compare(int natoms, int n1, rvec** eigvec1, int n2, rvec** eigvec2, real* eigval1, int neig1, real* eigval2, int neig2)
233 double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
235 n = std::min(n1, n2);
237 n = std::min(n, std::min(neig1, neig2));
238 fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
241 for (i = 0; i < n; i++)
248 eigval1[i] = std::sqrt(eigval1[i]);
251 for (i = n; i < neig1; i++)
253 trace1 += eigval1[i];
256 for (i = 0; i < n; i++)
263 eigval2[i] = std::sqrt(eigval2[i]);
268 // The static analyzer appears to be confused by the fact that the loop below
269 // starts from n instead of 0. However, given all the complex code it's
270 // better to be safe than sorry, so we check it with an assert.
271 // If we are in this comparison routine in the first place, neig2 should not be 0,
272 // so eigval2 should always be a valid pointer.
273 GMX_RELEASE_ASSERT(eigval2 != nullptr, "NULL pointer provided for eigval2");
275 for (i = n; i < neig2; i++)
277 trace2 += eigval2[i];
280 fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
281 if (neig1 != n || neig2 != n)
284 "this is %d%% and %d%% of the total trace\n",
285 gmx::roundToInt(100 * sum1 / trace1),
286 gmx::roundToInt(100 * sum2 / trace2));
288 fprintf(stdout, "Square root of the traces: %g and %g\n", std::sqrt(sum1), std::sqrt(sum2));
291 for (i = 0; i < n; i++)
294 for (j = 0; j < n; j++)
297 for (k = 0; k < natoms; k++)
299 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
301 tmp += eigval2[j] * ip * ip;
303 sab += eigval1[i] * tmp;
306 samsb2 = sum1 + sum2 - 2 * sab;
312 fprintf(stdout, "The overlap of the covariance matrices:\n");
313 fprintf(stdout, " normalized: %.3f\n", 1 - std::sqrt(samsb2 / (sum1 + sum2)));
314 tmp = 1 - sab / std::sqrt(sum1 * sum2);
319 fprintf(stdout, " shape: %.3f\n", 1 - std::sqrt(tmp));
323 static void inprod_matrix(const char* matfile,
337 int i, x1, y1, x, y, nlevels;
339 real inp, *t_x, *t_y, maxval;
347 for (y1 = 0; y1 < nx; y1++)
349 if (outvec[y1] < nvec2)
351 t_y[ny] = eignr2[outvec[y1]] + 1;
360 for (y = 0; y < ny; y++)
362 t_y[y] = eignr2[y] + 1;
366 fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n", nx, nvec2);
371 for (x1 = 0; x1 < nx; x1++)
382 t_x[x1] = eignr1[x] + 1;
383 fprintf(stderr, " %d", eignr1[x] + 1);
384 for (y1 = 0; y1 < ny; y1++)
389 while (outvec[y1] >= nvec2)
399 for (i = 0; i < natoms; i++)
401 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
403 mat[x1][y1] = std::abs(inp);
404 if (mat[x1][y1] > maxval)
406 maxval = mat[x1][y1];
410 fprintf(stderr, "\n");
418 out = gmx_ffopen(matfile, "w");
419 write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2", nx, ny, t_x, t_y, mat, 0.0, maxval, rlo, rhi, &nlevels);
423 static void overlap(const char* outfile,
431 const gmx_output_env_t* oenv)
437 fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
438 for (i = 0; i < noutvec; i++)
440 fprintf(stderr, "%d ", outvec[i] + 1);
442 fprintf(stderr, "\n");
444 out = xvgropen(outfile, "Subspace overlap", "Eigenvectors of trajectory 2", "Overlap", oenv);
445 if (output_env_get_print_xvgr_codes(oenv))
447 fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
450 for (x = 0; x < nvec2; x++)
452 for (v = 0; v < noutvec; v++)
456 for (i = 0; i < natoms; i++)
458 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
460 overlap += gmx::square(inp);
462 fprintf(out, "%5d %5.3f\n", eignr2[x] + 1, overlap / noutvec);
468 static void project(const char* trajfile,
469 const t_topology* top,
472 const char* projfile,
473 const char* twodplotfile,
474 const char* threedplotfile,
475 const char* filterfile,
477 const char* extremefile,
481 const t_atoms* atoms,
496 const gmx_output_env_t* oenv)
498 FILE* xvgrout = nullptr;
499 int nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
500 t_trxstatus* out = nullptr;
502 int noutvec_extr, imin, imax;
507 real t, inp, **inprod = nullptr;
508 char str[STRLEN], str2[STRLEN], *c;
511 gmx_rmpbc_t gpbc = nullptr;
517 noutvec_extr = noutvec;
527 snew(inprod, noutvec + 1);
531 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n", filterfile);
532 for (i = 0; i < noutvec; i++)
534 fprintf(stderr, "%d ", outvec[i] + 1);
536 fprintf(stderr, "\n");
537 out = open_trx(filterfile, "w");
542 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
546 "the number of atoms in your trajectory (%d) is larger than the number of "
547 "atoms in your structure file (%d)",
555 gpbc = gmx_rmpbc_init(&top->idef, pbcType, nat);
558 for (i = 0; i < nat; i++)
568 gmx_rmpbc(gpbc, nat, box, xread);
570 if (nframes >= snew_size)
573 for (i = 0; i < noutvec + 1; i++)
575 srenew(inprod[i], snew_size);
578 inprod[noutvec][nframes] = t;
579 /* calculate x: a fitted struture of the selected atoms */
582 reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
583 do_fit(nat, w_rls, xref, xread);
585 for (i = 0; i < natoms; i++)
587 copy_rvec(xread[index[i]], x[i]);
590 for (v = 0; v < noutvec; v++)
593 /* calculate (mass-weighted) projection */
595 for (i = 0; i < natoms; i++)
597 inp += (eigvec[vec][i][0] * (x[i][0] - xav[i][0])
598 + eigvec[vec][i][1] * (x[i][1] - xav[i][1])
599 + eigvec[vec][i][2] * (x[i][2] - xav[i][2]))
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, nullptr, nullptr);
624 } while (read_next_x(oenv, status, &t, xread, box));
634 snew(xread, atoms->nr);
639 gmx_rmpbc_done(gpbc);
645 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL if projfile is non-NULL");
646 snew(ylabel, noutvec);
647 for (v = 0; v < noutvec; v++)
649 sprintf(str, "vec %d", eignr[outvec[v]] + 1);
650 ylabel[v] = gmx_strdup(str);
652 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
653 write_xvgr_graphs(projfile,
658 output_env_get_xvgr_tlabel(oenv),
664 output_env_get_time_factor(oenv),
672 sprintf(str, "projection on eigenvector %d (%s)", eignr[outvec[0]] + 1, proj_unit);
673 sprintf(str2, "projection on eigenvector %d (%s)", eignr[outvec[noutvec - 1]] + 1, proj_unit);
674 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2, oenv);
675 for (i = 0; i < nframes; i++)
677 if (bSplit && i > 0 && std::abs(inprod[noutvec][i]) < 1e-5)
679 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
681 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec - 1][i]);
698 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
702 bPDB = fn2ftp(threedplotfile) == efPDB;
704 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
706 b4D = bPDB && (noutvec >= 4);
710 "You have selected four or more eigenvectors:\n"
711 "fourth eigenvector will be plotted "
712 "in bfactor field of pdb file\n");
714 "4D proj. of traj. on eigenv. %d, %d, %d and %d",
715 eignr[outvec[0]] + 1,
716 eignr[outvec[1]] + 1,
717 eignr[outvec[2]] + 1,
718 eignr[outvec[3]] + 1);
723 "3D proj. of traj. on eigenv. %d, %d and %d",
724 eignr[outvec[0]] + 1,
725 eignr[outvec[1]] + 1,
726 eignr[outvec[2]] + 1);
728 init_t_atoms(&atoms, nframes, FALSE);
731 atnm = gmx_strdup("C");
732 resnm = gmx_strdup("PRJ");
736 fact = 10000.0 / nframes;
743 for (i = 0; i < nframes; i++)
745 atoms.atomname[i] = &atnm;
746 atoms.atom[i].resind = i;
747 atoms.resinfo[i].name = &resnm;
748 atoms.resinfo[i].nr = static_cast<int>(std::ceil(i * fact));
749 atoms.resinfo[i].ic = ' ';
750 x[i][XX] = inprod[0][i];
751 x[i][YY] = inprod[1][i];
752 x[i][ZZ] = inprod[2][i];
758 if ((b4D || bSplit) && bPDB)
760 GMX_RELEASE_ASSERT(inprod != nullptr,
761 "inprod must be non-NULL with 4D or split PDB output options");
763 out = gmx_ffopen(threedplotfile, "w");
764 fprintf(out, "HEADER %s\n", str);
767 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
770 for (i = 0; i < atoms.nr; i++)
772 if (j > 0 && bSplit && std::abs(inprod[noutvec][i]) < 1e-5)
774 fprintf(out, "TER\n");
777 gmx_fprintf_pdb_atomline(out,
794 fprintf(out, "CONECT%5d%5d\n", i, i + 1);
798 fprintf(out, "TER\n");
803 write_sto_conf(threedplotfile, str, &atoms, x, nullptr, pbcType, box);
810 snew(pmin, noutvec_extr);
811 snew(pmax, noutvec_extr);
814 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL");
815 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
816 fprintf(stderr, "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
819 for (v = 0; v < noutvec_extr; v++)
821 for (i = 0; i < nframes; i++)
823 if (inprod[v][i] < inprod[v][imin])
827 if (inprod[v][i] > inprod[v][imax])
832 pmin[v] = inprod[v][imin];
833 pmax[v] = inprod[v][imax];
834 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n", eignr[outvec[v]] + 1, pmin[v], imin, pmax[v], imax);
842 /* build format string for filename: */
843 std::strcpy(str, extremefile); /* copy filename */
844 c = std::strrchr(str, '.'); /* find where extention begins */
845 std::strcpy(str2, c); /* get extention */
846 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
847 for (v = 0; v < noutvec_extr; v++)
849 /* make filename using format string */
850 if (noutvec_extr == 1)
852 std::strcpy(str2, extremefile);
856 sprintf(str2, str, eignr[outvec[v]] + 1);
858 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n", nextr, outvec[v] + 1, str2);
859 out = open_trx(str2, "w");
860 for (frame = 0; frame < nextr; frame++)
862 if ((extreme == 0) && (nextr <= 3))
864 for (i = 0; i < natoms; i++)
866 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
869 for (i = 0; i < natoms; i++)
871 for (d = 0; d < DIM; d++)
875 + (pmin[v] * (nextr - frame - 1) + pmax[v] * frame) / (nextr - 1)
876 * eigvec[outvec[v]][i][d] / sqrtm[i]);
879 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, nullptr, nullptr);
886 fprintf(stderr, "\n");
889 static void components(const char* outfile,
895 const gmx_output_env_t* oenv)
902 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
904 snew(ylabel, noutvec);
907 for (i = 0; i < natoms; i++)
911 for (g = 0; g < noutvec; g++)
914 sprintf(str, "vec %d", eignr[v] + 1);
915 ylabel[g] = gmx_strdup(str);
917 for (s = 0; s < 4; s++)
919 snew(y[g][s], natoms);
921 for (i = 0; i < natoms; i++)
923 y[g][0][i] = norm(eigvec[v][i]);
924 for (s = 0; s < 3; s++)
926 y[g][s + 1][i] = eigvec[v][i][s];
930 write_xvgr_graphs(outfile,
933 "Eigenvector components",
934 "black: total, red: x, green: y, blue: z",
945 fprintf(stderr, "\n");
948 static void rmsf(const char* outfile,
957 const gmx_output_env_t* oenv)
964 for (i = 0; i < neig; i++)
972 fprintf(stderr, "Writing rmsf to %s\n", outfile);
974 snew(ylabel, noutvec);
977 for (i = 0; i < natoms; i++)
981 for (g = 0; g < noutvec; g++)
984 if (eignr[v] >= neig)
987 "Selected vector %d is larger than the number of eigenvalues (%d)",
991 sprintf(str, "vec %d", eignr[v] + 1);
992 ylabel[g] = gmx_strdup(str);
994 for (i = 0; i < natoms; i++)
996 y[g][i] = std::sqrt(eigval[eignr[v]] * norm2(eigvec[v][i])) / sqrtm[i];
1000 outfile, noutvec, 1, "RMS fluctuation (nm) ", nullptr, "Atom number", ylabel, natoms, x, y, nullptr, 1, TRUE, FALSE, oenv);
1001 fprintf(stderr, "\n");
1004 int gmx_anaeig(int argc, char* argv[])
1006 static const char* desc[] = {
1007 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
1008 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
1009 "([gmx-nmeig]).[PAR]",
1011 "When a trajectory is projected on eigenvectors, all structures are",
1012 "fitted to the structure in the eigenvector file, if present, otherwise",
1013 "to the structure in the structure file. When no run input file is",
1014 "supplied, periodicity will not be taken into account. Most analyses",
1015 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
1016 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
1018 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
1019 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
1021 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
1022 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
1024 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
1025 "[TT]-first[tt] to [TT]-last[tt].",
1026 "The projections of a trajectory on the eigenvectors of its",
1027 "covariance matrix are called principal components (pc's).",
1028 "It is often useful to check the cosine content of the pc's,",
1029 "since the pc's of random diffusion are cosines with the number",
1030 "of periods equal to half the pc index.",
1031 "The cosine content of the pc's can be calculated with the program",
1032 "[gmx-analyze].[PAR]",
1034 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
1035 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
1037 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
1038 "three selected eigenvectors.[PAR]",
1040 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
1041 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
1043 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
1044 "on the average structure and interpolate [TT]-nframes[tt] frames",
1045 "between them, or set your own extremes with [TT]-max[tt]. The",
1046 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
1047 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1048 "will be written to separate files. Chain identifiers will be added",
1049 "when writing a [REF].pdb[ref] file with two or three structures (you",
1050 "can use [TT]rasmol -nmrpdb[tt] to view such a [REF].pdb[ref] file).[PAR]",
1052 "Overlap calculations between covariance analysis",
1053 "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^",
1055 "[BB]Note:[bb] the analysis should use the same fitting structure",
1057 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1058 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1059 "in file [TT]-v[tt].[PAR]",
1061 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1062 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1063 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1064 "have been set explicitly.[PAR]",
1066 "When [TT]-v[tt] and [TT]-v2[tt] are given, a single number for the",
1067 "overlap between the covariance matrices is generated. Note that the",
1068 "eigenvalues are by default read from the timestamp field in the",
1069 "eigenvector input files, but when [TT]-eig[tt], or [TT]-eig2[tt] are",
1070 "given, the corresponding eigenvalues are used instead. The formulas are::",
1072 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))",
1073 " normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))",
1074 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))",
1076 "where M1 and M2 are the two covariance matrices and tr is the trace",
1077 "of a matrix. The numbers are proportional to the overlap of the square",
1078 "root of the fluctuations. The normalized overlap is the most useful",
1079 "number, it is 1 for identical matrices and 0 when the sampled",
1080 "subspaces are orthogonal.[PAR]",
1081 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1082 "computed based on the Quasiharmonic approach and based on",
1083 "Schlitter's formula."
1085 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1086 static real max = 0.0, temp = 298.15;
1087 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1089 { "-first", FALSE, etINT, { &first }, "First eigenvector for analysis (-1 is select)" },
1090 { "-last", FALSE, etINT, { &last }, "Last eigenvector for analysis (-1 is till the last)" },
1091 { "-skip", FALSE, etINT, { &skip }, "Only analyse every nr-th frame" },
1096 "Maximum for projection of the eigenvector on the average structure, "
1097 "max=0 gives the extremes" },
1098 { "-nframes", FALSE, etINT, { &nextr }, "Number of frames for the extremes output" },
1099 { "-split", FALSE, etBOOL, { &bSplit }, "Split eigenvector projections where time is zero" },
1104 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1105 { "-temp", FALSE, etREAL, { &temp }, "Temperature for entropy calculations" },
1110 "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic "
1111 "approximation. When you do a rotational and/or translational fit prior to the "
1112 "covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which "
1113 "should not be taken into account when computing the entropy." }
1115 #define NPA asize(pa)
1118 PbcType pbcType = PbcType::Unset;
1119 const t_atoms* atoms = nullptr;
1120 rvec * xtop, *xref1, *xref2, *xrefp = nullptr;
1121 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1122 int nvec1, nvec2, *eignr1 = nullptr, *eignr2 = nullptr;
1123 rvec * xav1, *xav2, **eigvec1 = nullptr, **eigvec2 = nullptr;
1125 real totmass, *sqrtm, *w_rls, t;
1128 const char* indexfile;
1130 int nout, *iout, noutvec, *outvec, nfit;
1131 int * index = nullptr, *ifit = nullptr;
1132 const char * VecFile, *Vec2File, *topfile;
1133 const char * EigFile, *Eig2File;
1134 const char * CompFile, *RmsfFile, *ProjOnVecFile;
1135 const char * TwoDPlotFile, *ThreeDPlotFile;
1136 const char * FilterFile, *ExtremeFile;
1137 const char * OverlapFile, *InpMatFile;
1138 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1139 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1140 real * eigval1 = nullptr, *eigval2 = nullptr;
1143 gmx_output_env_t* oenv;
1147 { efTRN, "-v", "eigenvec", ffREAD }, { efTRN, "-v2", "eigenvec2", ffOPTRD },
1148 { efTRX, "-f", nullptr, ffOPTRD }, { efTPS, nullptr, nullptr, ffOPTRD },
1149 { efNDX, nullptr, nullptr, ffOPTRD }, { efXVG, "-eig", "eigenval", ffOPTRD },
1150 { efXVG, "-eig2", "eigenval2", ffOPTRD }, { efXVG, "-comp", "eigcomp", ffOPTWR },
1151 { efXVG, "-rmsf", "eigrmsf", ffOPTWR }, { efXVG, "-proj", "proj", ffOPTWR },
1152 { efXVG, "-2d", "2dproj", ffOPTWR }, { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1153 { efTRX, "-filt", "filtered", ffOPTWR }, { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1154 { efXVG, "-over", "overlap", ffOPTWR }, { efXPM, "-inpr", "inprod", ffOPTWR }
1156 #define NFILE asize(fnm)
1158 if (!parse_common_args(
1159 &argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW, NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1164 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1166 VecFile = opt2fn("-v", NFILE, fnm);
1167 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1168 topfile = ftp2fn(efTPS, NFILE, fnm);
1169 EigFile = opt2fn_null("-eig", NFILE, fnm);
1170 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1171 CompFile = opt2fn_null("-comp", NFILE, fnm);
1172 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1173 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1174 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1175 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1176 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1177 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1178 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1179 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1181 bProj = (ProjOnVecFile != nullptr) || (TwoDPlotFile != nullptr) || (ThreeDPlotFile != nullptr)
1182 || (FilterFile != nullptr) || (ExtremeFile != nullptr);
1183 bFirstLastSet = opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1184 bFirstToLast = (CompFile != nullptr) || (RmsfFile != nullptr) || (ProjOnVecFile != nullptr)
1185 || (FilterFile != nullptr) || (OverlapFile != nullptr)
1186 || (((ExtremeFile != nullptr) || (InpMatFile != nullptr)) && bFirstLastSet);
1187 bVec2 = (Vec2File != nullptr) || (OverlapFile != nullptr) || (InpMatFile != nullptr);
1188 bM = (RmsfFile != nullptr) || bProj;
1189 bTraj = (ProjOnVecFile != nullptr) || (FilterFile != nullptr)
1190 || ((ExtremeFile != nullptr) && (max == 0)) || (TwoDPlotFile != nullptr)
1191 || (ThreeDPlotFile != nullptr);
1192 bIndex = bM || bProj;
1193 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj || (FilterFile != nullptr)
1194 || (bIndex && (indexfile != nullptr));
1195 bCompare = (Vec2File != nullptr) || (Eig2File != nullptr);
1196 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1199 VecFile, &natoms, &bFit1, &xref1, &bDMR1, &xav1, &bDMA1, &nvec1, &eignr1, &eigvec1, &eigval1);
1200 neig1 = std::min(nvec1, DIM * natoms);
1201 if (nvec1 != DIM * natoms)
1204 "Warning: number of eigenvectors %d does not match three times\n"
1205 "the number of atoms %d in %s. Using %d eigenvectors.\n\n",
1212 /* Overwrite eigenvalues from separate files if the user provides them */
1213 if (EigFile != nullptr)
1215 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1216 if (neig_tmp != neig1)
1219 "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",
1224 srenew(eigval1, neig1);
1225 for (j = 0; j < neig1; j++)
1227 real tmp = eigval1[j];
1228 eigval1[j] = xvgdata[1][j];
1229 if (debug && (eigval1[j] != tmp))
1231 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n", j, tmp, eigval1[j]);
1234 for (j = 0; j < i; j++)
1239 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1247 "Can not calculate entropies from mass-weighted eigenvalues, redo the "
1248 "analysis without mass-weighting");
1250 printf("The Entropy due to the Schlitter formula is %g J/mol K\n",
1251 calcSchlitterEntropy(gmx::arrayRefFromArray(eigval1, neig1), temp, FALSE));
1252 printf("The Entropy due to the Quasiharmonic analysis is %g J/mol K\n",
1253 calcQuasiHarmonicEntropy(gmx::arrayRefFromArray(eigval1, neig1), temp, FALSE, 1.0));
1260 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1264 Vec2File, &natoms2, &bFit2, &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1266 neig2 = std::min(nvec2, DIM * natoms2);
1269 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1278 if (Eig2File != nullptr)
1280 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1281 srenew(eigval2, neig2);
1282 for (j = 0; j < neig2; j++)
1284 eigval2[j] = xvgdata[1][j];
1286 for (j = 0; j < i; j++)
1291 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1295 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1299 if ((xref1 == nullptr) && (bM || bTraj))
1315 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &pbcType, &xtop, nullptr, topbox, bM);
1317 gpbc = gmx_rmpbc_init(&top.idef, pbcType, atoms->nr);
1318 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1319 /* Fitting is only required for the projection */
1322 if (xref1 == nullptr)
1324 printf("\nNote: the structure in %s should be the same\n"
1325 " as the one used for the fit in g_covar\n",
1328 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1329 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1331 snew(w_rls, atoms->nr);
1332 for (i = 0; (i < nfit); i++)
1336 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1340 w_rls[ifit[i]] = 1.0;
1344 snew(xrefp, atoms->nr);
1345 if (xref1 != nullptr)
1347 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1351 "you selected a group with %d elements instead of %d, your selection "
1352 "does not fit the reference structure in the eigenvector file.",
1356 for (i = 0; (i < nfit); i++)
1358 copy_rvec(xref1[i], xrefp[ifit[i]]);
1363 /* The top coordinates are the fitting reference */
1364 for (i = 0; (i < nfit); i++)
1366 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1368 reset_x(nfit, ifit, atoms->nr, nullptr, xrefp, w_rls);
1371 gmx_rmpbc_done(gpbc);
1376 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1377 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1380 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1385 snew(sqrtm, natoms);
1388 proj_unit = "u\\S1/2\\Nnm";
1389 for (i = 0; (i < natoms); i++)
1391 sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
1397 for (i = 0; (i < natoms); i++)
1407 for (i = 0; (i < natoms); i++)
1409 for (d = 0; (d < DIM); d++)
1411 t += gmx::square((xav1[i][d] - xav2[i][d]) * sqrtm[i]);
1412 totmass += gmx::square(sqrtm[i]);
1416 "RMSD (without fit) between the two average structures:"
1418 std::sqrt(t / totmass));
1423 last = natoms * DIM;
1429 /* make an index from first to last */
1430 nout = last - first + 1;
1432 for (i = 0; i < nout; i++)
1434 iout[i] = first - 1 + i;
1437 else if (ThreeDPlotFile)
1439 /* make an index of first+(0,1,2) and last */
1440 nout = bPDB3D ? 4 : 3;
1441 nout = std::min(last - first + 1, nout);
1443 iout[0] = first - 1;
1447 iout[2] = first + 1;
1449 iout[nout - 1] = last - 1;
1453 /* make an index of first and last */
1456 iout[0] = first - 1;
1462 printf("Select eigenvectors for output, end your selection with 0\n");
1469 srenew(iout, nout + 1);
1470 if (1 != scanf("%d", &iout[nout]))
1472 gmx_fatal(FARGS, "Error reading user input");
1475 } while (iout[nout] >= 0);
1479 /* make an index of the eigenvectors which are present */
1482 for (i = 0; i < nout; i++)
1485 while ((j < nvec1) && (eignr1[j] != iout[i]))
1489 if ((j < nvec1) && (eignr1[j] == iout[i]))
1491 outvec[noutvec] = j;
1495 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1498 fprintf(stderr, ":");
1499 for (j = 0; j < noutvec; j++)
1501 fprintf(stderr, " %d", eignr1[outvec[j]] + 1);
1504 fprintf(stderr, "\n");
1508 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1513 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1, neig1, oenv);
1518 project(bTraj ? opt2fn("-f", NFILE, fnm) : nullptr,
1519 bTop ? &top : nullptr,
1551 overlap(OverlapFile, natoms, eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1557 InpMatFile, natoms, nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2, bFirstLastSet, noutvec, outvec);
1562 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1566 if (!CompFile && !bProj && !OverlapFile && !InpMatFile && !bCompare && !bEntropy)
1569 "\nIf you want some output,"
1570 " set one (or two or ...) of the output file options\n");
1574 view_all(oenv, NFILE, fnm);