src/gromacs/gmxana/gmx_anaeig.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <cassert>
  40 #include <cmath>
  41 #include <cstdlib>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45 #include <string>
  46
  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"
  69
  70 #include "thermochemistry.h"
  71
  72 static const char *proj_unit;
  73
  74 static real tick_spacing(real range, int minticks)
  75 {
  76     real sp;
  77
  78     if (range <= 0)
  79     {
  80         return 1.0;
  81     }
  82
  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)
  85     {
  86         sp = sp/2;
  87     }
  88
  89     return sp;
  90 }
  91
  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)
  98 {
  99     FILE *out;
 100     int   g, s, i;
 101     real  ymin, ymax, xsp, ysp;
 102
 103     out = gmx_ffopen(file, "w");
 104     if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
 105     {
 106         fprintf(out, "@ autoscale onread none\n");
 107     }
 108     for (g = 0; g < ngraphs; g++)
 109     {
 110         if (y)
 111         {
 112             ymin = y[g][0];
 113             ymax = y[g][0];
 114             for (i = 0; i < n; i++)
 115             {
 116                 if (y[g][i] < ymin)
 117                 {
 118                     ymin = y[g][i];
 119                 }
 120                 if (y[g][i] > ymax)
 121                 {
 122                     ymax = y[g][i];
 123                 }
 124             }
 125         }
 126         else
 127         {
 128             assert(sy);
 129             ymin = sy[g][0][0];
 130             ymax = sy[g][0][0];
 131             for (s = 0; s < nsetspergraph; s++)
 132             {
 133                 for (i = 0; i < n; i++)
 134                 {
 135                     if (sy[g][s][i] < ymin)
 136                     {
 137                         ymin = sy[g][s][i];
 138                     }
 139                     if (sy[g][s][i] > ymax)
 140                     {
 141                         ymax = sy[g][s][i];
 142                     }
 143                 }
 144             }
 145         }
 146         if (bZero)
 147         {
 148             ymin = 0;
 149         }
 150         else
 151         {
 152             ymin = ymin-0.1*(ymax-ymin);
 153         }
 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))
 158         {
 159             fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
 160             if (g == 0)
 161             {
 162                 fprintf(out, "@ title \"%s\"\n", title);
 163                 if (subtitle)
 164                 {
 165                     fprintf(out, "@ subtitle \"%s\"\n", subtitle);
 166                 }
 167             }
 168             if (g == ngraphs-1)
 169             {
 170                 fprintf(out, "@ xaxis  label \"%s\"\n", xlabel.c_str());
 171             }
 172             else
 173             {
 174                 fprintf(out, "@ xaxis  ticklabel off\n");
 175             }
 176             if (n > 1)
 177             {
 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);
 182             }
 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))
 197             {
 198                 fprintf(out, "@ zeroxaxis bar on\n");
 199                 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
 200             }
 201         }
 202         for (s = 0; s < nsetspergraph; s++)
 203         {
 204             for (i = 0; i < n; i++)
 205             {
 206                 if (bSplit && i > 0 && std::abs(x[i]) < 1e-5)
 207                 {
 208                     fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 209                 }
 210                 fprintf(out, "%10.4f %10.5f\n",
 211                         x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
 212             }
 213             fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 214         }
 215     }
 216     gmx_ffclose(out);
 217 }
 218
 219 static void
 220 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
 221         real *eigval1, int neig1, real *eigval2, int neig2)
 222 {
 223     int    n;
 224     int    i, j, k;
 225     double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
 226
 227     n = std::min(n1, n2);
 228
 229     n = std::min(n, std::min(neig1, neig2));
 230     fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
 231
 232     sum1 = 0;
 233     for (i = 0; i < n; i++)
 234     {
 235         if (eigval1[i] < 0)
 236         {
 237             eigval1[i] = 0;
 238         }
 239         sum1      += eigval1[i];
 240         eigval1[i] = std::sqrt(eigval1[i]);
 241     }
 242     trace1 = sum1;
 243     for (i = n; i < neig1; i++)
 244     {
 245         trace1 += eigval1[i];
 246     }
 247     sum2 = 0;
 248     for (i = 0; i < n; i++)
 249     {
 250         if (eigval2[i] < 0)
 251         {
 252             eigval2[i] = 0;
 253         }
 254         sum2      += eigval2[i];
 255         eigval2[i] = std::sqrt(eigval2[i]);
 256     }
 257     trace2 = sum2;
 258
 259
 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");
 266
 267     for (i = n; i < neig2; i++)
 268     {
 269         trace2 += eigval2[i];
 270     }
 271
 272     fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
 273     if (neig1 != n || neig2 != n)
 274     {
 275         fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
 276                 gmx::roundToInt(100*sum1/trace1), gmx::roundToInt(100*sum2/trace2));
 277     }
 278     fprintf(stdout, "Square root of the traces: %g and %g\n",
 279             std::sqrt(sum1), std::sqrt(sum2));
 280
 281     sab = 0;
 282     for (i = 0; i < n; i++)
 283     {
 284         tmp = 0;
 285         for (j = 0; j < n; j++)
 286         {
 287             ip = 0;
 288             for (k = 0; k < natoms; k++)
 289             {
 290                 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
 291             }
 292             tmp += eigval2[j]*ip*ip;
 293         }
 294         sab += eigval1[i]*tmp;
 295     }
 296
 297     samsb2 = sum1+sum2-2*sab;
 298     if (samsb2 < 0)
 299     {
 300         samsb2 = 0;
 301     }
 302
 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)
 307     {
 308         tmp = 0;
 309     }
 310     fprintf(stdout, "       shape:  %.3f\n", 1-std::sqrt(tmp));
 311 }
 312
 313
 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)
 318 {
 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;
 325
 326     snew(t_y, nvec2);
 327     if (bSelect)
 328     {
 329         nx = noutvec;
 330         ny = 0;
 331         for (y1 = 0; y1 < nx; y1++)
 332         {
 333             if (outvec[y1] < nvec2)
 334             {
 335                 t_y[ny] = eignr2[outvec[y1]]+1;
 336                 ny++;
 337             }
 338         }
 339     }
 340     else
 341     {
 342         nx = nvec1;
 343         ny = nvec2;
 344         for (y = 0; y < ny; y++)
 345         {
 346             t_y[y] = eignr2[y]+1;
 347         }
 348     }
 349
 350     fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
 351             nx, nvec2);
 352
 353     snew(mat, nx);
 354     snew(t_x, nx);
 355     maxval = 0;
 356     for (x1 = 0; x1 < nx; x1++)
 357     {
 358         snew(mat[x1], ny);
 359         if (bSelect)
 360         {
 361             x = outvec[x1];
 362         }
 363         else
 364         {
 365             x = x1;
 366         }
 367         t_x[x1] = eignr1[x]+1;
 368         fprintf(stderr, " %d", eignr1[x]+1);
 369         for (y1 = 0; y1 < ny; y1++)
 370         {
 371             inp = 0;
 372             if (bSelect)
 373             {
 374                 while (outvec[y1] >= nvec2)
 375                 {
 376                     y1++;
 377                 }
 378                 y = outvec[y1];
 379             }
 380             else
 381             {
 382                 y = y1;
 383             }
 384             for (i = 0; i < natoms; i++)
 385             {
 386                 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
 387             }
 388             mat[x1][y1] = std::abs(inp);
 389             if (mat[x1][y1] > maxval)
 390             {
 391                 maxval = mat[x1][y1];
 392             }
 393         }
 394     }
 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);
 403 }
 404
 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)
 410 {
 411     FILE *out;
 412     int   i, v, vec, x;
 413     real  overlap, inp;
 414
 415     fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
 416     for (i = 0; i < noutvec; i++)
 417     {
 418         fprintf(stderr, "%d ", outvec[i]+1);
 419     }
 420     fprintf(stderr, "\n");
 421
 422     out = xvgropen(outfile, "Subspace overlap",
 423                    "Eigenvectors of trajectory 2", "Overlap", oenv);
 424     if (output_env_get_print_xvgr_codes(oenv))
 425     {
 426         fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
 427     }
 428     overlap = 0;
 429     for (x = 0; x < nvec2; x++)
 430     {
 431         for (v = 0; v < noutvec; v++)
 432         {
 433             vec = outvec[v];
 434             inp = 0;
 435             for (i = 0; i < natoms; i++)
 436             {
 437                 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
 438             }
 439             overlap += gmx::square(inp);
 440         }
 441         fprintf(out, "%5d  %5.3f\n", eignr2[x]+1, overlap/noutvec);
 442     }
 443
 444     xvgrclose(out);
 445 }
 446
 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)
 457 {
 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], *c;
 469     const char **ylabel;
 470     real         fact;
 471     gmx_rmpbc_t  gpbc = nullptr;
 472
 473     snew(x, natoms);
 474
 475     if (bExtrAll)
 476     {
 477         noutvec_extr = noutvec;
 478     }
 479     else
 480     {
 481         noutvec_extr = 1;
 482     }
 483
 484
 485     if (trajfile)
 486     {
 487         snew(inprod, noutvec+1);
 488
 489         if (filterfile)
 490         {
 491             fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
 492                     filterfile);
 493             for (i = 0; i < noutvec; i++)
 494             {
 495                 fprintf(stderr, "%d ", outvec[i]+1);
 496             }
 497             fprintf(stderr, "\n");
 498             out = open_trx(filterfile, "w");
 499         }
 500         snew_size = 0;
 501         nfr       = 0;
 502         nframes   = 0;
 503         nat       = read_first_x(oenv, &status, trajfile, &t, &xread, box);
 504         if (nat > atoms->nr)
 505         {
 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);
 507         }
 508         snew(all_at, nat);
 509
 510         if (top)
 511         {
 512             gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
 513         }
 514
 515         for (i = 0; i < nat; i++)
 516         {
 517             all_at[i] = i;
 518         }
 519         do
 520         {
 521             if (nfr % skip == 0)
 522             {
 523                 if (top)
 524                 {
 525                     gmx_rmpbc(gpbc, nat, box, xread);
 526                 }
 527                 if (nframes >= snew_size)
 528                 {
 529                     snew_size += 100;
 530                     for (i = 0; i < noutvec+1; i++)
 531                     {
 532                         srenew(inprod[i], snew_size);
 533                     }
 534                 }
 535                 inprod[noutvec][nframes] = t;
 536                 /* calculate x: a fitted struture of the selected atoms */
 537                 if (bFit)
 538                 {
 539                     reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
 540                     do_fit(nat, w_rls, xref, xread);
 541                 }
 542                 for (i = 0; i < natoms; i++)
 543                 {
 544                     copy_rvec(xread[index[i]], x[i]);
 545                 }
 546
 547                 for (v = 0; v < noutvec; v++)
 548                 {
 549                     vec = outvec[v];
 550                     /* calculate (mass-weighted) projection */
 551                     inp = 0;
 552                     for (i = 0; i < natoms; i++)
 553                     {
 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];
 557                     }
 558                     inprod[v][nframes] = inp;
 559                 }
 560                 if (filterfile)
 561                 {
 562                     for (i = 0; i < natoms; i++)
 563                     {
 564                         for (d = 0; d < DIM; d++)
 565                         {
 566                             /* misuse xread for output */
 567                             xread[index[i]][d] = xav[i][d];
 568                             for (v = 0; v < noutvec; v++)
 569                             {
 570                                 xread[index[i]][d] +=
 571                                     inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
 572                             }
 573                         }
 574                     }
 575                     write_trx(out, natoms, index, atoms, 0, t, box, xread, nullptr, nullptr);
 576                 }
 577                 nframes++;
 578             }
 579             nfr++;
 580         }
 581         while (read_next_x(oenv, status, &t, xread, box));
 582         close_trx(status);
 583         sfree(x);
 584         if (filterfile)
 585         {
 586             close_trx(out);
 587         }
 588     }
 589     else
 590     {
 591         snew(xread, atoms->nr);
 592     }
 593
 594     if (top)
 595     {
 596         gmx_rmpbc_done(gpbc);
 597     }
 598
 599
 600     if (projfile)
 601     {
 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++)
 605         {
 606             sprintf(str, "vec %d", eignr[outvec[v]]+1);
 607             ylabel[v] = gmx_strdup(str);
 608         }
 609         sprintf(str, "projection on eigenvectors (%s)", proj_unit);
 610         write_xvgr_graphs(projfile, noutvec, 1, str, nullptr, output_env_get_xvgr_tlabel(oenv),
 611                           ylabel,
 612                           nframes, inprod[noutvec], inprod, nullptr,
 613                           output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
 614     }
 615
 616     if (twodplotfile)
 617     {
 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,
 623                            oenv);
 624         for (i = 0; i < nframes; i++)
 625         {
 626             if (bSplit && i > 0 && std::abs(inprod[noutvec][i]) < 1e-5)
 627             {
 628                 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 629             }
 630             fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
 631         }
 632         xvgrclose(xvgrout);
 633     }
 634
 635     if (threedplotfile)
 636     {
 637         t_atoms     atoms;
 638         rvec       *x;
 639         real       *b = nullptr;
 640         matrix      box;
 641         char       *resnm, *atnm;
 642         gmx_bool    bPDB, b4D;
 643         FILE       *out;
 644
 645         if (noutvec < 3)
 646         {
 647             gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
 648         }
 649
 650         /* initialize */
 651         bPDB = fn2ftp(threedplotfile) == efPDB;
 652         clear_mat(box);
 653         box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
 654
 655         b4D = bPDB && (noutvec >= 4);
 656         if (b4D)
 657         {
 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);
 664         }
 665         else
 666         {
 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);
 669         }
 670         init_t_atoms(&atoms, nframes, FALSE);
 671         snew(x, nframes);
 672         snew(b, nframes);
 673         atnm  = gmx_strdup("C");
 674         resnm = gmx_strdup("PRJ");
 675
 676         if (nframes > 10000)
 677         {
 678             fact = 10000.0/nframes;
 679         }
 680         else
 681         {
 682             fact = 1.0;
 683         }
 684
 685         for (i = 0; i < nframes; i++)
 686         {
 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];
 695             if (b4D)
 696             {
 697                 b[i]  = inprod[3][i];
 698             }
 699         }
 700         if ( ( b4D || bSplit ) && bPDB)
 701         {
 702             GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL with 4D or split PDB output options");
 703
 704             out = gmx_ffopen(threedplotfile, "w");
 705             fprintf(out, "HEADER    %s\n", str);
 706             if (b4D)
 707             {
 708                 fprintf(out, "REMARK    %s\n", "fourth dimension plotted as B-factor");
 709             }
 710             j = 0;
 711             for (i = 0; i < atoms.nr; i++)
 712             {
 713                 if (j > 0 && bSplit && std::abs(inprod[noutvec][i]) < 1e-5)
 714                 {
 715                     fprintf(out, "TER\n");
 716                     j = 0;
 717                 }
 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], "");
 720                 if (j > 0)
 721                 {
 722                     fprintf(out, "CONECT%5d%5d\n", i, i+1);
 723                 }
 724                 j++;
 725             }
 726             fprintf(out, "TER\n");
 727             gmx_ffclose(out);
 728         }
 729         else
 730         {
 731             write_sto_conf(threedplotfile, str, &atoms, x, nullptr, ePBC, box);
 732         }
 733         done_atom(&atoms);
 734     }
 735
 736     if (extremefile)
 737     {
 738         snew(pmin, noutvec_extr);
 739         snew(pmax, noutvec_extr);
 740         if (extreme == 0)
 741         {
 742             GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL");
 743             fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
 744             fprintf(stderr,
 745                     "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
 746             imin = 0;
 747             imax = 0;
 748             for (v = 0; v < noutvec_extr; v++)
 749             {
 750                 for (i = 0; i < nframes; i++)
 751                 {
 752                     if (inprod[v][i] < inprod[v][imin])
 753                     {
 754                         imin = i;
 755                     }
 756                     if (inprod[v][i] > inprod[v][imax])
 757                     {
 758                         imax = i;
 759                     }
 760                 }
 761                 pmin[v] = inprod[v][imin];
 762                 pmax[v] = inprod[v][imax];
 763                 fprintf(stderr, "%7d     %10.6f %10d %10.6f %10d\n",
 764                         eignr[outvec[v]]+1,
 765                         pmin[v], imin, pmax[v], imax);
 766             }
 767         }
 768         else
 769         {
 770             pmin[0] = -extreme;
 771             pmax[0] =  extreme;
 772         }
 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++)
 779         {
 780             /* make filename using format string */
 781             if (noutvec_extr == 1)
 782             {
 783                 std::strcpy(str2, extremefile);
 784             }
 785             else
 786             {
 787                 sprintf(str2, str, eignr[outvec[v]]+1);
 788             }
 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++)
 793             {
 794                 if ((extreme == 0) && (nextr <= 3))
 795                 {
 796                     for (i = 0; i < natoms; i++)
 797                     {
 798                         atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
 799                     }
 800                 }
 801                 for (i = 0; i < natoms; i++)
 802                 {
 803                     for (d = 0; d < DIM; d++)
 804                     {
 805                         xread[index[i]][d] =
 806                             (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
 807                              *eigvec[outvec[v]][i][d]/sqrtm[i]);
 808                     }
 809                 }
 810                 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, nullptr, nullptr);
 811             }
 812             close_trx(out);
 813         }
 814         sfree(pmin);
 815         sfree(pmax);
 816     }
 817     fprintf(stderr, "\n");
 818 }
 819
 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)
 824 {
 825     int         g, s, v, i;
 826     real       *x, ***y;
 827     char        str[STRLEN];
 828     const char**ylabel;
 829
 830     fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
 831
 832     snew(ylabel, noutvec);
 833     snew(y, noutvec);
 834     snew(x, natoms);
 835     for (i = 0; i < natoms; i++)
 836     {
 837         x[i] = i+1;
 838     }
 839     for (g = 0; g < noutvec; g++)
 840     {
 841         v = outvec[g];
 842         sprintf(str, "vec %d", eignr[v]+1);
 843         ylabel[g] = gmx_strdup(str);
 844         snew(y[g], 4);
 845         for (s = 0; s < 4; s++)
 846         {
 847             snew(y[g][s], natoms);
 848         }
 849         for (i = 0; i < natoms; i++)
 850         {
 851             y[g][0][i] = norm(eigvec[v][i]);
 852             for (s = 0; s < 3; s++)
 853             {
 854                 y[g][s+1][i] = eigvec[v][i][s];
 855             }
 856         }
 857     }
 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");
 863 }
 864
 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)
 870 {
 871     int          g, v, i;
 872     real        *x, **y;
 873     char         str[STRLEN];
 874     const char **ylabel;
 875
 876     for (i = 0; i < neig; i++)
 877     {
 878         if (eigval[i] < 0)
 879         {
 880             eigval[i] = 0;
 881         }
 882     }
 883
 884     fprintf(stderr, "Writing rmsf to %s\n", outfile);
 885
 886     snew(ylabel, noutvec);
 887     snew(y, noutvec);
 888     snew(x, natoms);
 889     for (i = 0; i < natoms; i++)
 890     {
 891         x[i] = i+1;
 892     }
 893     for (g = 0; g < noutvec; g++)
 894     {
 895         v = outvec[g];
 896         if (eignr[v] >= neig)
 897         {
 898             gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
 899         }
 900         sprintf(str, "vec %d", eignr[v]+1);
 901         ylabel[g] = gmx_strdup(str);
 902         snew(y[g], natoms);
 903         for (i = 0; i < natoms; i++)
 904         {
 905             y[g][i] = std::sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
 906         }
 907     }
 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");
 912 }
 913
 914 int gmx_anaeig(int argc, char *argv[])
 915 {
 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]",
 920
 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]",
 927
 928         "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
 929         "[TT]-first[tt] to [TT]-last[tt].[PAR]",
 930
 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]",
 933
 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]",
 943
 944         "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
 945         "[TT]-first[tt] and [TT]-last[tt].[PAR]",
 946
 947         "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
 948         "three selected eigenvectors.[PAR]",
 949
 950         "[TT]-filt[tt]: filter the trajectory to show only the motion along",
 951         "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
 952
 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]",
 961
 962         "Overlap calculations between covariance analysis",
 963         "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^",
 964         "",
 965         "[BB]Note:[bb] the analysis should use the same fitting structure",
 966         "",
 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]",
 970
 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]",
 975
 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::",
 981         "",
 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))",
 985         "",
 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."
 994     };
 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;
 998     t_pargs            pa[]   = {
 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." }
1018     };
1019 #define NPA asize(pa)
1020
1021     t_topology        top;
1022     int               ePBC  = -1;
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;
1028     matrix            topbox;
1029     real              totmass, *sqrtm, *w_rls, t;
1030     int               natoms;
1031     char             *grpname;
1032     const char       *indexfile;
1033     int               i, j, d;
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;
1045     int               neig1, neig2;
1046     double          **xvgdata;
1047     gmx_output_env_t *oenv;
1048     gmx_rmpbc_t       gpbc;
1049
1050     t_filenm          fnm[] = {
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 }
1067     };
1068 #define NFILE asize(fnm)
1069
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))
1073     {
1074         return 0;
1075     }
1076
1077     indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1078
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);
1093
1094     bProj  = (ProjOnVecFile != nullptr) || (TwoDPlotFile != nullptr) || (ThreeDPlotFile != nullptr)
1095         || (FilterFile != nullptr) || (ExtremeFile != nullptr);
1096     bFirstLastSet  =
1097         opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1098     bFirstToLast = (CompFile != nullptr) || (RmsfFile != nullptr) || (ProjOnVecFile != nullptr) || (FilterFile != nullptr) ||
1099         (OverlapFile != nullptr) || (((ExtremeFile != nullptr) || (InpMatFile != nullptr)) && bFirstLastSet);
1100     bVec2  = (Vec2File != nullptr) || (OverlapFile != nullptr) || (InpMatFile != nullptr);
1101     bM     = (RmsfFile != nullptr) || bProj;
1102     bTraj  = (ProjOnVecFile != nullptr) || (FilterFile != nullptr) || ((ExtremeFile != nullptr) && (max == 0))
1103         || (TwoDPlotFile != nullptr) || (ThreeDPlotFile != nullptr);
1104     bIndex = bM || bProj;
1105     bTPS   = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1106         (FilterFile != nullptr)  || (bIndex && (indexfile != nullptr));
1107     bCompare = (Vec2File != nullptr) || (Eig2File != nullptr);
1108     bPDB3D   = fn2ftp(ThreeDPlotFile) == efPDB;
1109
1110     read_eigenvectors(VecFile, &natoms, &bFit1,
1111                       &xref1, &bDMR1, &xav1, &bDMA1,
1112                       &nvec1, &eignr1, &eigvec1, &eigval1);
1113     neig1 = std::min(nvec1, DIM*natoms);
1114     if (nvec1 != DIM*natoms)
1115     {
1116         fprintf(stderr, "Warning: number of eigenvectors %d does not match three times\n"
1117                 "the number of atoms %d in %s. Using %d eigenvectors.\n\n",
1118                 nvec1, natoms, VecFile, neig1);
1119     }
1120
1121     /* Overwrite eigenvalues from separate files if the user provides them */
1122     if (EigFile != nullptr)
1123     {
1124         int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1125         if (neig_tmp != neig1)
1126         {
1127             fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1128         }
1129         neig1 = neig_tmp;
1130         srenew(eigval1, neig1);
1131         for (j = 0; j < neig1; j++)
1132         {
1133             real tmp = eigval1[j];
1134             eigval1[j] = xvgdata[1][j];
1135             if (debug && (eigval1[j] != tmp))
1136             {
1137                 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1138                         j, tmp, eigval1[j]);
1139             }
1140         }
1141         for (j = 0; j < i; j++)
1142         {
1143             sfree(xvgdata[j]);
1144         }
1145         sfree(xvgdata);
1146         fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1147     }
1148
1149     if (bEntropy)
1150     {
1151         if (bDMA1)
1152         {
1153             gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1154         }
1155         printf("The Entropy due to the Schlitter formula is %g J/mol K\n",
1156                calcSchlitterEntropy(gmx::arrayRefFromArray(eigval1, neig1),
1157                                     temp, FALSE));
1158         printf("The Entropy due to the Quasiharmonic analysis is %g J/mol K\n",
1159                calcQuasiHarmonicEntropy(gmx::arrayRefFromArray(eigval1, neig1),
1160                                         temp, FALSE, 1.0));
1161     }
1162
1163     if (bVec2)
1164     {
1165         if (!Vec2File)
1166         {
1167             gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1168         }
1169         read_eigenvectors(Vec2File, &neig2, &bFit2,
1170                           &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1171
1172         neig2 = DIM*neig2;
1173         if (neig2 != neig1)
1174         {
1175             gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1176         }
1177     }
1178     else
1179     {
1180         nvec2 = 0;
1181         neig2 = 0;
1182     }
1183
1184     if (Eig2File != nullptr)
1185     {
1186         neig2 = read_xvg(Eig2File, &xvgdata, &i);
1187         srenew(eigval2, neig2);
1188         for (j = 0; j < neig2; j++)
1189         {
1190             eigval2[j] = xvgdata[1][j];
1191         }
1192         for (j = 0; j < i; j++)
1193         {
1194             sfree(xvgdata[j]);
1195         }
1196         sfree(xvgdata);
1197         fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1198     }
1199
1200
1201     if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1202     {
1203         bM = FALSE;
1204     }
1205     if ((xref1 == nullptr) && (bM || bTraj))
1206     {
1207         bTPS = TRUE;
1208     }
1209
1210     xtop  = nullptr;
1211     nfit  = 0;
1212     ifit  = nullptr;
1213     w_rls = nullptr;
1214
1215     if (!bTPS)
1216     {
1217         bTop = FALSE;
1218     }
1219     else
1220     {
1221         bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1222                              &top, &ePBC, &xtop, nullptr, topbox, bM);
1223         atoms = &top.atoms;
1224         gpbc  = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1225         gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1226         /* Fitting is only required for the projection */
1227         if (bProj && bFit1)
1228         {
1229             if (xref1 == nullptr)
1230             {
1231                 printf("\nNote: the structure in %s should be the same\n"
1232                        "      as the one used for the fit in g_covar\n", topfile);
1233             }
1234             printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1235             get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1236
1237             snew(w_rls, atoms->nr);
1238             for (i = 0; (i < nfit); i++)
1239             {
1240                 if (bDMR1)
1241                 {
1242                     w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1243                 }
1244                 else
1245                 {
1246                     w_rls[ifit[i]] = 1.0;
1247                 }
1248             }
1249
1250             snew(xrefp, atoms->nr);
1251             if (xref1 != nullptr)
1252             {
1253                 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1254                 if (natoms != nfit)
1255                 {
1256                     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);
1257                 }
1258                 for (i = 0; (i < nfit); i++)
1259                 {
1260                     copy_rvec(xref1[i], xrefp[ifit[i]]);
1261                 }
1262             }
1263             else
1264             {
1265                 /* The top coordinates are the fitting reference */
1266                 for (i = 0; (i < nfit); i++)
1267                 {
1268                     copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1269                 }
1270                 reset_x(nfit, ifit, atoms->nr, nullptr, xrefp, w_rls);
1271             }
1272         }
1273         gmx_rmpbc_done(gpbc);
1274     }
1275
1276     if (bIndex)
1277     {
1278         printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1279         get_index(atoms, indexfile, 1, &i, &index, &grpname);
1280         if (i != natoms)
1281         {
1282             gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1283         }
1284         printf("\n");
1285     }
1286
1287     snew(sqrtm, natoms);
1288     if (bM && bDMA1)
1289     {
1290         proj_unit = "u\\S1/2\\Nnm";
1291         for (i = 0; (i < natoms); i++)
1292         {
1293             sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
1294         }
1295     }
1296     else
1297     {
1298         proj_unit = "nm";
1299         for (i = 0; (i < natoms); i++)
1300         {
1301             sqrtm[i] = 1.0;
1302         }
1303     }
1304
1305     if (bVec2)
1306     {
1307         t       = 0;
1308         totmass = 0;
1309         for (i = 0; (i < natoms); i++)
1310         {
1311             for (d = 0; (d < DIM); d++)
1312             {
1313                 t       += gmx::square((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1314                 totmass += gmx::square(sqrtm[i]);
1315             }
1316         }
1317         fprintf(stdout, "RMSD (without fit) between the two average structures:"
1318                 " %.3f (nm)\n\n", std::sqrt(t/totmass));
1319     }
1320
1321     if (last == -1)
1322     {
1323         last = natoms*DIM;
1324     }
1325     if (first > -1)
1326     {
1327         if (bFirstToLast)
1328         {
1329             /* make an index from first to last */
1330             nout = last-first+1;
1331             snew(iout, nout);
1332             for (i = 0; i < nout; i++)
1333             {
1334                 iout[i] = first-1+i;
1335             }
1336         }
1337         else if (ThreeDPlotFile)
1338         {
1339             /* make an index of first+(0,1,2) and last */
1340             nout = bPDB3D ? 4 : 3;
1341             nout = std::min(last-first+1, nout);
1342             snew(iout, nout);
1343             iout[0] = first-1;
1344             iout[1] = first;
1345             if (nout > 3)
1346             {
1347                 iout[2] = first+1;
1348             }
1349             iout[nout-1] = last-1;
1350         }
1351         else
1352         {
1353             /* make an index of first and last */
1354             nout = 2;
1355             snew(iout, nout);
1356             iout[0] = first-1;
1357             iout[1] = last-1;
1358         }
1359     }
1360     else
1361     {
1362         printf("Select eigenvectors for output, end your selection with 0\n");
1363         nout = -1;
1364         iout = nullptr;
1365
1366         do
1367         {
1368             nout++;
1369             srenew(iout, nout+1);
1370             if (1 != scanf("%d", &iout[nout]))
1371             {
1372                 gmx_fatal(FARGS, "Error reading user input");
1373             }
1374             iout[nout]--;
1375         }
1376         while (iout[nout] >= 0);
1377
1378         printf("\n");
1379     }
1380     /* make an index of the eigenvectors which are present */
1381     snew(outvec, nout);
1382     noutvec = 0;
1383     for (i = 0; i < nout; i++)
1384     {
1385         j = 0;
1386         while ((j < nvec1) && (eignr1[j] != iout[i]))
1387         {
1388             j++;
1389         }
1390         if ((j < nvec1) && (eignr1[j] == iout[i]))
1391         {
1392             outvec[noutvec] = j;
1393             noutvec++;
1394         }
1395     }
1396     fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1397     if (noutvec <= 100)
1398     {
1399         fprintf(stderr, ":");
1400         for (j = 0; j < noutvec; j++)
1401         {
1402             fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1403         }
1404     }
1405     fprintf(stderr, "\n");
1406
1407     if (CompFile)
1408     {
1409         components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1410     }
1411
1412     if (RmsfFile)
1413     {
1414         rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1415              neig1, oenv);
1416     }
1417
1418     if (bProj)
1419     {
1420         project(bTraj ? opt2fn("-f", NFILE, fnm) : nullptr,
1421                 bTop ? &top : nullptr, ePBC, topbox,
1422                 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1423                 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1424                 bFit1, xrefp, nfit, ifit, w_rls,
1425                 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1426                 oenv);
1427     }
1428
1429     if (OverlapFile)
1430     {
1431         overlap(OverlapFile, natoms,
1432                 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1433     }
1434
1435     if (InpMatFile)
1436     {
1437         inprod_matrix(InpMatFile, natoms,
1438                       nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1439                       bFirstLastSet, noutvec, outvec);
1440     }
1441
1442     if (bCompare)
1443     {
1444         compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1445     }
1446
1447
1448     if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1449         !bCompare && !bEntropy)
1450     {
1451         fprintf(stderr, "\nIf you want some output,"
1452                 " set one (or two or ...) of the output file options\n");
1453     }
1454
1455
1456     view_all(oenv, NFILE, fnm);
1457
1458     return 0;
1459 }