Renamed entropy.* to thermochemistry.*
[gromacs.git] / src / gromacs / gmxana / gmx_anaeig.cpp
blob3366b3a6930ce8b7939b92a9e5639f2d548cfa0a
1 /*
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.
37 #include "gmxpre.h"
39 #include <cassert>
40 #include <cmath>
41 #include <cstdlib>
42 #include <cstring>
44 #include <algorithm>
45 #include <string>
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 const char *proj_unit;
74 static real tick_spacing(real range, int minticks)
76 real sp;
78 if (range <= 0)
80 return 1.0;
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)
86 sp = sp/2;
89 return sp;
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)
99 FILE *out;
100 int g, s, i;
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++)
110 if (y)
112 ymin = y[g][0];
113 ymax = y[g][0];
114 for (i = 0; i < n; i++)
116 if (y[g][i] < ymin)
118 ymin = y[g][i];
120 if (y[g][i] > ymax)
122 ymax = y[g][i];
126 else
128 assert(sy);
129 ymin = sy[g][0][0];
130 ymax = sy[g][0][0];
131 for (s = 0; s < nsetspergraph; s++)
133 for (i = 0; i < n; i++)
135 if (sy[g][s][i] < ymin)
137 ymin = sy[g][s][i];
139 if (sy[g][s][i] > ymax)
141 ymax = sy[g][s][i];
146 if (bZero)
148 ymin = 0;
150 else
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);
160 if (g == 0)
162 fprintf(out, "@ title \"%s\"\n", title);
163 if (subtitle)
165 fprintf(out, "@ subtitle \"%s\"\n", subtitle);
168 if (g == ngraphs-1)
170 fprintf(out, "@ xaxis label \"%s\"\n", xlabel.c_str());
172 else
174 fprintf(out, "@ xaxis ticklabel off\n");
176 if (n > 1)
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) ? "&" : "");
216 gmx_ffclose(out);
219 static void
220 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
221 real *eigval1, int neig1, real *eigval2, int neig2)
223 int n;
224 int i, j, k;
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);
232 sum1 = 0;
233 for (i = 0; i < n; i++)
235 if (eigval1[i] < 0)
237 eigval1[i] = 0;
239 sum1 += eigval1[i];
240 eigval1[i] = std::sqrt(eigval1[i]);
242 trace1 = sum1;
243 for (i = n; i < neig1; i++)
245 trace1 += eigval1[i];
247 sum2 = 0;
248 for (i = 0; i < n; i++)
250 if (eigval2[i] < 0)
252 eigval2[i] = 0;
254 sum2 += eigval2[i];
255 eigval2[i] = std::sqrt(eigval2[i]);
257 trace2 = sum2;
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));
281 sab = 0;
282 for (i = 0; i < n; i++)
284 tmp = 0;
285 for (j = 0; j < n; j++)
287 ip = 0;
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;
298 if (samsb2 < 0)
300 samsb2 = 0;
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);
306 if (tmp < 0)
308 tmp = 0;
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, int *eignr2, rvec **eigvec2,
317 gmx_bool bSelect, int noutvec, int *outvec)
319 FILE *out;
320 real **mat;
321 int i, x1, y1, x, y, nlevels;
322 int nx, ny;
323 real inp, *t_x, *t_y, maxval;
324 t_rgb rlo, rhi;
326 snew(t_y, nvec2);
327 if (bSelect)
329 nx = noutvec;
330 ny = 0;
331 for (y1 = 0; y1 < nx; y1++)
333 if (outvec[y1] < nvec2)
335 t_y[ny] = eignr2[outvec[y1]]+1;
336 ny++;
340 else
342 nx = nvec1;
343 ny = nvec2;
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",
351 nx, nvec2);
353 snew(mat, nx);
354 snew(t_x, nx);
355 maxval = 0;
356 for (x1 = 0; x1 < nx; x1++)
358 snew(mat[x1], ny);
359 if (bSelect)
361 x = outvec[x1];
363 else
365 x = x1;
367 t_x[x1] = eignr1[x]+1;
368 fprintf(stderr, " %d", eignr1[x]+1);
369 for (y1 = 0; y1 < ny; y1++)
371 inp = 0;
372 if (bSelect)
374 while (outvec[y1] >= nvec2)
376 y1++;
378 y = outvec[y1];
380 else
382 y = y1;
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;
398 nlevels = 41;
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);
402 gmx_ffclose(out);
405 static void overlap(const char *outfile, int natoms,
406 rvec **eigvec1,
407 int nvec2, int *eignr2, rvec **eigvec2,
408 int noutvec, int *outvec,
409 const gmx_output_env_t *oenv)
411 FILE *out;
412 int i, v, vec, x;
413 real overlap, inp;
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);
428 overlap = 0;
429 for (x = 0; x < nvec2; x++)
431 for (v = 0; v < noutvec; v++)
433 vec = outvec[v];
434 inp = 0;
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);
444 xvgrclose(out);
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 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;
461 t_trxstatus *status;
462 int noutvec_extr, imin, imax;
463 real *pmin, *pmax;
464 int *all_at;
465 matrix box;
466 rvec *xread, *x;
467 real t, inp, **inprod = nullptr;
468 char str[STRLEN], str2[STRLEN], **ylabel, *c;
469 real fact;
470 gmx_rmpbc_t gpbc = nullptr;
472 snew(x, natoms);
474 if (bExtrAll)
476 noutvec_extr = noutvec;
478 else
480 noutvec_extr = 1;
484 if (trajfile)
486 snew(inprod, noutvec+1);
488 if (filterfile)
490 fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
491 filterfile);
492 for (i = 0; i < noutvec; i++)
494 fprintf(stderr, "%d ", outvec[i]+1);
496 fprintf(stderr, "\n");
497 out = open_trx(filterfile, "w");
499 snew_size = 0;
500 nfr = 0;
501 nframes = 0;
502 nat = read_first_x(oenv, &status, trajfile, &t, &xread, box);
503 if (nat > atoms->nr)
505 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);
507 snew(all_at, nat);
509 if (top)
511 gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
514 for (i = 0; i < nat; i++)
516 all_at[i] = i;
520 if (nfr % skip == 0)
522 if (top)
524 gmx_rmpbc(gpbc, nat, box, xread);
526 if (nframes >= snew_size)
528 snew_size += 100;
529 for (i = 0; i < noutvec+1; i++)
531 srenew(inprod[i], snew_size);
534 inprod[noutvec][nframes] = t;
535 /* calculate x: a fitted struture of the selected atoms */
536 if (bFit)
538 reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
539 do_fit(nat, w_rls, xref, xread);
541 for (i = 0; i < natoms; i++)
543 copy_rvec(xread[index[i]], x[i]);
546 for (v = 0; v < noutvec; v++)
548 vec = outvec[v];
549 /* calculate (mass-weighted) projection */
550 inp = 0;
551 for (i = 0; i < natoms; i++)
553 inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
554 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
555 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
557 inprod[v][nframes] = inp;
559 if (filterfile)
561 for (i = 0; i < natoms; i++)
563 for (d = 0; d < DIM; d++)
565 /* misuse xread for output */
566 xread[index[i]][d] = xav[i][d];
567 for (v = 0; v < noutvec; v++)
569 xread[index[i]][d] +=
570 inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
574 write_trx(out, natoms, index, atoms, 0, t, box, xread, nullptr, nullptr);
576 nframes++;
578 nfr++;
580 while (read_next_x(oenv, status, &t, xread, box));
581 close_trx(status);
582 sfree(x);
583 if (filterfile)
585 close_trx(out);
588 else
590 snew(xread, atoms->nr);
593 if (top)
595 gmx_rmpbc_done(gpbc);
599 if (projfile)
601 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL if projfile is non-NULL");
602 snew(ylabel, noutvec);
603 for (v = 0; v < noutvec; v++)
605 sprintf(str, "vec %d", eignr[outvec[v]]+1);
606 ylabel[v] = gmx_strdup(str);
608 sprintf(str, "projection on eigenvectors (%s)", proj_unit);
609 write_xvgr_graphs(projfile, noutvec, 1, str, nullptr, output_env_get_xvgr_tlabel(oenv),
610 (const char **)ylabel,
611 nframes, inprod[noutvec], inprod, nullptr,
612 output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
615 if (twodplotfile)
617 sprintf(str, "projection on eigenvector %d (%s)",
618 eignr[outvec[0]]+1, proj_unit);
619 sprintf(str2, "projection on eigenvector %d (%s)",
620 eignr[outvec[noutvec-1]]+1, proj_unit);
621 xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
622 oenv);
623 for (i = 0; i < nframes; i++)
625 if (bSplit && i > 0 && std::abs(inprod[noutvec][i]) < 1e-5)
627 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
629 fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
631 xvgrclose(xvgrout);
634 if (threedplotfile)
636 t_atoms atoms;
637 rvec *x;
638 real *b = nullptr;
639 matrix box;
640 char *resnm, *atnm;
641 gmx_bool bPDB, b4D;
642 FILE *out;
644 if (noutvec < 3)
646 gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
649 /* initialize */
650 bPDB = fn2ftp(threedplotfile) == efPDB;
651 clear_mat(box);
652 box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
654 b4D = bPDB && (noutvec >= 4);
655 if (b4D)
657 fprintf(stderr, "You have selected four or more eigenvectors:\n"
658 "fourth eigenvector will be plotted "
659 "in bfactor field of pdb file\n");
660 sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
661 eignr[outvec[0]]+1, eignr[outvec[1]]+1,
662 eignr[outvec[2]]+1, eignr[outvec[3]]+1);
664 else
666 sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
667 eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
669 init_t_atoms(&atoms, nframes, FALSE);
670 snew(x, nframes);
671 snew(b, nframes);
672 atnm = gmx_strdup("C");
673 resnm = gmx_strdup("PRJ");
675 if (nframes > 10000)
677 fact = 10000.0/nframes;
679 else
681 fact = 1.0;
684 for (i = 0; i < nframes; i++)
686 atoms.atomname[i] = &atnm;
687 atoms.atom[i].resind = i;
688 atoms.resinfo[i].name = &resnm;
689 atoms.resinfo[i].nr = static_cast<int>(std::ceil(i*fact));
690 atoms.resinfo[i].ic = ' ';
691 x[i][XX] = inprod[0][i];
692 x[i][YY] = inprod[1][i];
693 x[i][ZZ] = inprod[2][i];
694 if (b4D)
696 b[i] = inprod[3][i];
699 if ( ( b4D || bSplit ) && bPDB)
701 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL with 4D or split PDB output options");
703 out = gmx_ffopen(threedplotfile, "w");
704 fprintf(out, "HEADER %s\n", str);
705 if (b4D)
707 fprintf(out, "REMARK %s\n", "fourth dimension plotted as B-factor");
709 j = 0;
710 for (i = 0; i < atoms.nr; i++)
712 if (j > 0 && bSplit && std::abs(inprod[noutvec][i]) < 1e-5)
714 fprintf(out, "TER\n");
715 j = 0;
717 gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
718 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
719 if (j > 0)
721 fprintf(out, "CONECT%5d%5d\n", i, i+1);
723 j++;
725 fprintf(out, "TER\n");
726 gmx_ffclose(out);
728 else
730 write_sto_conf(threedplotfile, str, &atoms, x, nullptr, ePBC, box);
732 done_atom(&atoms);
735 if (extremefile)
737 snew(pmin, noutvec_extr);
738 snew(pmax, noutvec_extr);
739 if (extreme == 0)
741 GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL");
742 fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
743 fprintf(stderr,
744 "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
745 imin = 0;
746 imax = 0;
747 for (v = 0; v < noutvec_extr; v++)
749 for (i = 0; i < nframes; i++)
751 if (inprod[v][i] < inprod[v][imin])
753 imin = i;
755 if (inprod[v][i] > inprod[v][imax])
757 imax = i;
760 pmin[v] = inprod[v][imin];
761 pmax[v] = inprod[v][imax];
762 fprintf(stderr, "%7d %10.6f %10d %10.6f %10d\n",
763 eignr[outvec[v]]+1,
764 pmin[v], imin, pmax[v], imax);
767 else
769 pmin[0] = -extreme;
770 pmax[0] = extreme;
772 /* build format string for filename: */
773 std::strcpy(str, extremefile); /* copy filename */
774 c = std::strrchr(str, '.'); /* find where extention begins */
775 std::strcpy(str2, c); /* get extention */
776 sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
777 for (v = 0; v < noutvec_extr; v++)
779 /* make filename using format string */
780 if (noutvec_extr == 1)
782 std::strcpy(str2, extremefile);
784 else
786 sprintf(str2, str, eignr[outvec[v]]+1);
788 fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
789 nextr, outvec[v]+1, str2);
790 out = open_trx(str2, "w");
791 for (frame = 0; frame < nextr; frame++)
793 if ((extreme == 0) && (nextr <= 3))
795 for (i = 0; i < natoms; i++)
797 atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
800 for (i = 0; i < natoms; i++)
802 for (d = 0; d < DIM; d++)
804 xread[index[i]][d] =
805 (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
806 *eigvec[outvec[v]][i][d]/sqrtm[i]);
809 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, nullptr, nullptr);
811 close_trx(out);
813 sfree(pmin);
814 sfree(pmax);
816 fprintf(stderr, "\n");
819 static void components(const char *outfile, int natoms,
820 int *eignr, rvec **eigvec,
821 int noutvec, int *outvec,
822 const gmx_output_env_t *oenv)
824 int g, s, v, i;
825 real *x, ***y;
826 char str[STRLEN], **ylabel;
828 fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
830 snew(ylabel, noutvec);
831 snew(y, noutvec);
832 snew(x, natoms);
833 for (i = 0; i < natoms; i++)
835 x[i] = i+1;
837 for (g = 0; g < noutvec; g++)
839 v = outvec[g];
840 sprintf(str, "vec %d", eignr[v]+1);
841 ylabel[g] = gmx_strdup(str);
842 snew(y[g], 4);
843 for (s = 0; s < 4; s++)
845 snew(y[g][s], natoms);
847 for (i = 0; i < natoms; i++)
849 y[g][0][i] = norm(eigvec[v][i]);
850 for (s = 0; s < 3; s++)
852 y[g][s+1][i] = eigvec[v][i][s];
856 write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
857 "black: total, red: x, green: y, blue: z",
858 "Atom number", (const char **)ylabel,
859 natoms, x, nullptr, y, 1, FALSE, FALSE, oenv);
860 fprintf(stderr, "\n");
863 static void rmsf(const char *outfile, int natoms, real *sqrtm,
864 int *eignr, rvec **eigvec,
865 int noutvec, int *outvec,
866 real *eigval, int neig,
867 const gmx_output_env_t *oenv)
869 int g, v, i;
870 real *x, **y;
871 char str[STRLEN], **ylabel;
873 for (i = 0; i < neig; i++)
875 if (eigval[i] < 0)
877 eigval[i] = 0;
881 fprintf(stderr, "Writing rmsf to %s\n", outfile);
883 snew(ylabel, noutvec);
884 snew(y, noutvec);
885 snew(x, natoms);
886 for (i = 0; i < natoms; i++)
888 x[i] = i+1;
890 for (g = 0; g < noutvec; g++)
892 v = outvec[g];
893 if (eignr[v] >= neig)
895 gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
897 sprintf(str, "vec %d", eignr[v]+1);
898 ylabel[g] = gmx_strdup(str);
899 snew(y[g], natoms);
900 for (i = 0; i < natoms; i++)
902 y[g][i] = std::sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
905 write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", nullptr,
906 "Atom number", (const char **)ylabel,
907 natoms, x, y, nullptr, 1, TRUE, FALSE, oenv);
908 fprintf(stderr, "\n");
911 int gmx_anaeig(int argc, char *argv[])
913 static const char *desc[] = {
914 "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
915 "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
916 "([gmx-nmeig]).[PAR]",
918 "When a trajectory is projected on eigenvectors, all structures are",
919 "fitted to the structure in the eigenvector file, if present, otherwise",
920 "to the structure in the structure file. When no run input file is",
921 "supplied, periodicity will not be taken into account. Most analyses",
922 "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
923 "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
925 "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
926 "[TT]-first[tt] to [TT]-last[tt].[PAR]",
928 "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
929 "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
931 "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
932 "[TT]-first[tt] to [TT]-last[tt].",
933 "The projections of a trajectory on the eigenvectors of its",
934 "covariance matrix are called principal components (pc's).",
935 "It is often useful to check the cosine content of the pc's,",
936 "since the pc's of random diffusion are cosines with the number",
937 "of periods equal to half the pc index.",
938 "The cosine content of the pc's can be calculated with the program",
939 "[gmx-analyze].[PAR]",
941 "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
942 "[TT]-first[tt] and [TT]-last[tt].[PAR]",
944 "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
945 "three selected eigenvectors.[PAR]",
947 "[TT]-filt[tt]: filter the trajectory to show only the motion along",
948 "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
950 "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
951 "on the average structure and interpolate [TT]-nframes[tt] frames",
952 "between them, or set your own extremes with [TT]-max[tt]. The",
953 "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
954 "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
955 "will be written to separate files. Chain identifiers will be added",
956 "when writing a [REF].pdb[ref] file with two or three structures (you",
957 "can use [TT]rasmol -nmrpdb[tt] to view such a [REF].pdb[ref] file).[PAR]",
959 "Overlap calculations between covariance analysis",
960 "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^",
962 "[BB]Note:[bb] the analysis should use the same fitting structure",
964 "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
965 "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
966 "in file [TT]-v[tt].[PAR]",
968 "[TT]-inpr[tt]: calculate a matrix of inner-products between",
969 "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
970 "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
971 "have been set explicitly.[PAR]",
973 "When [TT]-v[tt] and [TT]-v2[tt] are given, a single number for the",
974 "overlap between the covariance matrices is generated. Note that the",
975 "eigenvalues are by default read from the timestamp field in the",
976 "eigenvector input files, but when [TT]-eig[tt], or [TT]-eig2[tt] are",
977 "given, the corresponding eigenvalues are used instead. The formulas are::",
979 " difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))",
980 " normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))",
981 " shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))",
983 "where M1 and M2 are the two covariance matrices and tr is the trace",
984 "of a matrix. The numbers are proportional to the overlap of the square",
985 "root of the fluctuations. The normalized overlap is the most useful",
986 "number, it is 1 for identical matrices and 0 when the sampled",
987 "subspaces are orthogonal.[PAR]",
988 "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
989 "computed based on the Quasiharmonic approach and based on",
990 "Schlitter's formula."
992 static int first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
993 static real max = 0.0, temp = 298.15;
994 static gmx_bool bSplit = FALSE, bEntropy = FALSE;
995 t_pargs pa[] = {
996 { "-first", FALSE, etINT, {&first},
997 "First eigenvector for analysis (-1 is select)" },
998 { "-last", FALSE, etINT, {&last},
999 "Last eigenvector for analysis (-1 is till the last)" },
1000 { "-skip", FALSE, etINT, {&skip},
1001 "Only analyse every nr-th frame" },
1002 { "-max", FALSE, etREAL, {&max},
1003 "Maximum for projection of the eigenvector on the average structure, "
1004 "max=0 gives the extremes" },
1005 { "-nframes", FALSE, etINT, {&nextr},
1006 "Number of frames for the extremes output" },
1007 { "-split", FALSE, etBOOL, {&bSplit},
1008 "Split eigenvector projections where time is zero" },
1009 { "-entropy", FALSE, etBOOL, {&bEntropy},
1010 "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1011 { "-temp", FALSE, etREAL, {&temp},
1012 "Temperature for entropy calculations" },
1013 { "-nevskip", FALSE, etINT, {&nskip},
1014 "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." }
1016 #define NPA asize(pa)
1018 t_topology top;
1019 int ePBC = -1;
1020 const t_atoms *atoms = nullptr;
1021 rvec *xtop, *xref1, *xref2, *xrefp = nullptr;
1022 gmx_bool bDMR1, bDMA1, bDMR2, bDMA2;
1023 int nvec1, nvec2, *eignr1 = nullptr, *eignr2 = nullptr;
1024 rvec *xav1, *xav2, **eigvec1 = nullptr, **eigvec2 = nullptr;
1025 matrix topbox;
1026 real totmass, *sqrtm, *w_rls, t;
1027 int natoms;
1028 char *grpname;
1029 const char *indexfile;
1030 int i, j, d;
1031 int nout, *iout, noutvec, *outvec, nfit;
1032 int *index = nullptr, *ifit = nullptr;
1033 const char *VecFile, *Vec2File, *topfile;
1034 const char *EigFile, *Eig2File;
1035 const char *CompFile, *RmsfFile, *ProjOnVecFile;
1036 const char *TwoDPlotFile, *ThreeDPlotFile;
1037 const char *FilterFile, *ExtremeFile;
1038 const char *OverlapFile, *InpMatFile;
1039 gmx_bool bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1040 gmx_bool bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1041 real *eigval1 = nullptr, *eigval2 = nullptr;
1042 int neig1, neig2;
1043 double **xvgdata;
1044 gmx_output_env_t *oenv;
1045 gmx_rmpbc_t gpbc;
1047 t_filenm fnm[] = {
1048 { efTRN, "-v", "eigenvec", ffREAD },
1049 { efTRN, "-v2", "eigenvec2", ffOPTRD },
1050 { efTRX, "-f", nullptr, ffOPTRD },
1051 { efTPS, nullptr, nullptr, ffOPTRD },
1052 { efNDX, nullptr, nullptr, ffOPTRD },
1053 { efXVG, "-eig", "eigenval", ffOPTRD },
1054 { efXVG, "-eig2", "eigenval2", ffOPTRD },
1055 { efXVG, "-comp", "eigcomp", ffOPTWR },
1056 { efXVG, "-rmsf", "eigrmsf", ffOPTWR },
1057 { efXVG, "-proj", "proj", ffOPTWR },
1058 { efXVG, "-2d", "2dproj", ffOPTWR },
1059 { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1060 { efTRX, "-filt", "filtered", ffOPTWR },
1061 { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1062 { efXVG, "-over", "overlap", ffOPTWR },
1063 { efXPM, "-inpr", "inprod", ffOPTWR }
1065 #define NFILE asize(fnm)
1067 if (!parse_common_args(&argc, argv,
1068 PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
1069 NFILE, fnm, NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1071 return 0;
1074 indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1076 VecFile = opt2fn("-v", NFILE, fnm);
1077 Vec2File = opt2fn_null("-v2", NFILE, fnm);
1078 topfile = ftp2fn(efTPS, NFILE, fnm);
1079 EigFile = opt2fn_null("-eig", NFILE, fnm);
1080 Eig2File = opt2fn_null("-eig2", NFILE, fnm);
1081 CompFile = opt2fn_null("-comp", NFILE, fnm);
1082 RmsfFile = opt2fn_null("-rmsf", NFILE, fnm);
1083 ProjOnVecFile = opt2fn_null("-proj", NFILE, fnm);
1084 TwoDPlotFile = opt2fn_null("-2d", NFILE, fnm);
1085 ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1086 FilterFile = opt2fn_null("-filt", NFILE, fnm);
1087 ExtremeFile = opt2fn_null("-extr", NFILE, fnm);
1088 OverlapFile = opt2fn_null("-over", NFILE, fnm);
1089 InpMatFile = ftp2fn_null(efXPM, NFILE, fnm);
1091 bProj = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1092 || FilterFile || ExtremeFile;
1093 bFirstLastSet =
1094 opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1095 bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1096 OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1097 bVec2 = Vec2File || OverlapFile || InpMatFile;
1098 bM = RmsfFile || bProj;
1099 bTraj = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1100 || TwoDPlotFile || ThreeDPlotFile;
1101 bIndex = bM || bProj;
1102 bTPS = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1103 FilterFile || (bIndex && indexfile);
1104 bCompare = Vec2File || Eig2File;
1105 bPDB3D = fn2ftp(ThreeDPlotFile) == efPDB;
1107 read_eigenvectors(VecFile, &natoms, &bFit1,
1108 &xref1, &bDMR1, &xav1, &bDMA1,
1109 &nvec1, &eignr1, &eigvec1, &eigval1);
1110 neig1 = DIM*natoms;
1112 /* Overwrite eigenvalues from separate files if the user provides them */
1113 if (EigFile != nullptr)
1115 int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1116 if (neig_tmp != neig1)
1118 fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1120 neig1 = neig_tmp;
1121 srenew(eigval1, neig1);
1122 for (j = 0; j < neig1; j++)
1124 real tmp = eigval1[j];
1125 eigval1[j] = xvgdata[1][j];
1126 if (debug && (eigval1[j] != tmp))
1128 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1129 j, tmp, eigval1[j]);
1132 for (j = 0; j < i; j++)
1134 sfree(xvgdata[j]);
1136 sfree(xvgdata);
1137 fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1140 if (bEntropy)
1142 if (bDMA1)
1144 gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1146 printf("The Entropy due to the Schlitter formula is %g J/mol K\n",
1147 calc_entropy_schlitter(neig1, eigval1, temp, FALSE));
1150 if (bVec2)
1152 if (!Vec2File)
1154 gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1156 read_eigenvectors(Vec2File, &neig2, &bFit2,
1157 &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1159 neig2 = DIM*neig2;
1160 if (neig2 != neig1)
1162 gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1165 else
1167 nvec2 = 0;
1168 neig2 = 0;
1171 if (Eig2File != nullptr)
1173 neig2 = read_xvg(Eig2File, &xvgdata, &i);
1174 srenew(eigval2, neig2);
1175 for (j = 0; j < neig2; j++)
1177 eigval2[j] = xvgdata[1][j];
1179 for (j = 0; j < i; j++)
1181 sfree(xvgdata[j]);
1183 sfree(xvgdata);
1184 fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1188 if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1190 bM = FALSE;
1192 if ((xref1 == nullptr) && (bM || bTraj))
1194 bTPS = TRUE;
1197 xtop = nullptr;
1198 nfit = 0;
1199 ifit = nullptr;
1200 w_rls = nullptr;
1202 if (!bTPS)
1204 bTop = FALSE;
1206 else
1208 bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1209 &top, &ePBC, &xtop, nullptr, topbox, bM);
1210 atoms = &top.atoms;
1211 gpbc = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1212 gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1213 /* Fitting is only required for the projection */
1214 if (bProj && bFit1)
1216 if (xref1 == nullptr)
1218 printf("\nNote: the structure in %s should be the same\n"
1219 " as the one used for the fit in g_covar\n", topfile);
1221 printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1222 get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1224 snew(w_rls, atoms->nr);
1225 for (i = 0; (i < nfit); i++)
1227 if (bDMR1)
1229 w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1231 else
1233 w_rls[ifit[i]] = 1.0;
1237 snew(xrefp, atoms->nr);
1238 if (xref1 != nullptr)
1240 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1241 if (natoms != nfit)
1243 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);
1245 for (i = 0; (i < nfit); i++)
1247 copy_rvec(xref1[i], xrefp[ifit[i]]);
1250 else
1252 /* The top coordinates are the fitting reference */
1253 for (i = 0; (i < nfit); i++)
1255 copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1257 reset_x(nfit, ifit, atoms->nr, nullptr, xrefp, w_rls);
1260 gmx_rmpbc_done(gpbc);
1263 if (bIndex)
1265 printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1266 get_index(atoms, indexfile, 1, &i, &index, &grpname);
1267 if (i != natoms)
1269 gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1271 printf("\n");
1274 snew(sqrtm, natoms);
1275 if (bM && bDMA1)
1277 proj_unit = "u\\S1/2\\Nnm";
1278 for (i = 0; (i < natoms); i++)
1280 sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
1283 else
1285 proj_unit = "nm";
1286 for (i = 0; (i < natoms); i++)
1288 sqrtm[i] = 1.0;
1292 if (bVec2)
1294 t = 0;
1295 totmass = 0;
1296 for (i = 0; (i < natoms); i++)
1298 for (d = 0; (d < DIM); d++)
1300 t += gmx::square((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1301 totmass += gmx::square(sqrtm[i]);
1304 fprintf(stdout, "RMSD (without fit) between the two average structures:"
1305 " %.3f (nm)\n\n", std::sqrt(t/totmass));
1308 if (last == -1)
1310 last = natoms*DIM;
1312 if (first > -1)
1314 if (bFirstToLast)
1316 /* make an index from first to last */
1317 nout = last-first+1;
1318 snew(iout, nout);
1319 for (i = 0; i < nout; i++)
1321 iout[i] = first-1+i;
1324 else if (ThreeDPlotFile)
1326 /* make an index of first+(0,1,2) and last */
1327 nout = bPDB3D ? 4 : 3;
1328 nout = std::min(last-first+1, nout);
1329 snew(iout, nout);
1330 iout[0] = first-1;
1331 iout[1] = first;
1332 if (nout > 3)
1334 iout[2] = first+1;
1336 iout[nout-1] = last-1;
1338 else
1340 /* make an index of first and last */
1341 nout = 2;
1342 snew(iout, nout);
1343 iout[0] = first-1;
1344 iout[1] = last-1;
1347 else
1349 printf("Select eigenvectors for output, end your selection with 0\n");
1350 nout = -1;
1351 iout = nullptr;
1355 nout++;
1356 srenew(iout, nout+1);
1357 if (1 != scanf("%d", &iout[nout]))
1359 gmx_fatal(FARGS, "Error reading user input");
1361 iout[nout]--;
1363 while (iout[nout] >= 0);
1365 printf("\n");
1367 /* make an index of the eigenvectors which are present */
1368 snew(outvec, nout);
1369 noutvec = 0;
1370 for (i = 0; i < nout; i++)
1372 j = 0;
1373 while ((j < nvec1) && (eignr1[j] != iout[i]))
1375 j++;
1377 if ((j < nvec1) && (eignr1[j] == iout[i]))
1379 outvec[noutvec] = j;
1380 noutvec++;
1383 fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1384 if (noutvec <= 100)
1386 fprintf(stderr, ":");
1387 for (j = 0; j < noutvec; j++)
1389 fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1392 fprintf(stderr, "\n");
1394 if (CompFile)
1396 components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1399 if (RmsfFile)
1401 rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1402 neig1, oenv);
1405 if (bProj)
1407 project(bTraj ? opt2fn("-f", NFILE, fnm) : nullptr,
1408 bTop ? &top : nullptr, ePBC, topbox,
1409 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1410 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1411 bFit1, xrefp, nfit, ifit, w_rls,
1412 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1413 oenv);
1416 if (OverlapFile)
1418 overlap(OverlapFile, natoms,
1419 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1422 if (InpMatFile)
1424 inprod_matrix(InpMatFile, natoms,
1425 nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1426 bFirstLastSet, noutvec, outvec);
1429 if (bCompare)
1431 compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1435 if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1436 !bCompare && !bEntropy)
1438 fprintf(stderr, "\nIf you want some output,"
1439 " set one (or two or ...) of the output file options\n");
1443 view_all(oenv, NFILE, fnm);
1445 return 0;