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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include <cassert>
  41 #include <cmath>
  42 #include <cstdlib>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46 #include <string>
  47
  48 #include "gromacs/commandline/pargs.h"
  49 #include "gromacs/commandline/viewit.h"
  50 #include "gromacs/fileio/confio.h"
  51 #include "gromacs/fileio/matio.h"
  52 #include "gromacs/fileio/pdbio.h"
  53 #include "gromacs/fileio/tpxio.h"
  54 #include "gromacs/fileio/trxio.h"
  55 #include "gromacs/fileio/xvgr.h"
  56 #include "gromacs/gmxana/eigio.h"
  57 #include "gromacs/gmxana/gmx_ana.h"
  58 #include "gromacs/math/do_fit.h"
  59 #include "gromacs/math/functions.h"
  60 #include "gromacs/math/vec.h"
  61 #include "gromacs/pbcutil/rmpbc.h"
  62 #include "gromacs/topology/index.h"
  63 #include "gromacs/topology/topology.h"
  64 #include "gromacs/utility/arraysize.h"
  65 #include "gromacs/utility/cstringutil.h"
  66 #include "gromacs/utility/fatalerror.h"
  67 #include "gromacs/utility/futil.h"
  68 #include "gromacs/utility/gmxassert.h"
  69 #include "gromacs/utility/smalloc.h"
  70
  71 #include "thermochemistry.h"
  72
  73 static const char* proj_unit;
  74
  75 static real tick_spacing(real range, int minticks)
  76 {
  77     real sp;
  78
  79     if (range <= 0)
  80     {
  81         return 1.0;
  82     }
  83
  84     sp = 0.2 * std::exp(std::log(10.0) * std::ceil(std::log(range) / std::log(10.0)));
  85     while (range / sp < minticks - 1)
  86     {
  87         sp = sp / 2;
  88     }
  89
  90     return sp;
  91 }
  92
  93 static void write_xvgr_graphs(const char*             file,
  94                               int                     ngraphs,
  95                               int                     nsetspergraph,
  96                               const char*             title,
  97                               const char*             subtitle,
  98                               const std::string&      xlabel,
  99                               const char**            ylabel,
 100                               int                     n,
 101                               real*                   x,
 102                               real**                  y,
 103                               real***                 sy,
 104                               real                    scale_x,
 105                               gmx_bool                bZero,
 106                               gmx_bool                bSplit,
 107                               const gmx_output_env_t* oenv)
 108 {
 109     FILE* out;
 110     int   g, s, i;
 111     real  ymin, ymax, xsp, ysp;
 112
 113     out = gmx_ffopen(file, "w");
 114     if (output_env_get_xvg_format(oenv) == exvgXMGRACE)
 115     {
 116         fprintf(out, "@ autoscale onread none\n");
 117     }
 118     for (g = 0; g < ngraphs; g++)
 119     {
 120         if (y)
 121         {
 122             ymin = y[g][0];
 123             ymax = y[g][0];
 124             for (i = 0; i < n; i++)
 125             {
 126                 if (y[g][i] < ymin)
 127                 {
 128                     ymin = y[g][i];
 129                 }
 130                 if (y[g][i] > ymax)
 131                 {
 132                     ymax = y[g][i];
 133                 }
 134             }
 135         }
 136         else
 137         {
 138             assert(sy);
 139             ymin = sy[g][0][0];
 140             ymax = sy[g][0][0];
 141             for (s = 0; s < nsetspergraph; s++)
 142             {
 143                 for (i = 0; i < n; i++)
 144                 {
 145                     if (sy[g][s][i] < ymin)
 146                     {
 147                         ymin = sy[g][s][i];
 148                     }
 149                     if (sy[g][s][i] > ymax)
 150                     {
 151                         ymax = sy[g][s][i];
 152                     }
 153                 }
 154             }
 155         }
 156         if (bZero)
 157         {
 158             ymin = 0;
 159         }
 160         else
 161         {
 162             ymin = ymin - 0.1 * (ymax - ymin);
 163         }
 164         ymax = ymax + 0.1 * (ymax - ymin);
 165         xsp  = tick_spacing((x[n - 1] - x[0]) * scale_x, 4);
 166         ysp  = tick_spacing(ymax - ymin, 3);
 167         if (output_env_get_print_xvgr_codes(oenv))
 168         {
 169             fprintf(out, "@ with g%d\n@ g%d on\n", g, g);
 170             if (g == 0)
 171             {
 172                 fprintf(out, "@ title \"%s\"\n", title);
 173                 if (subtitle)
 174                 {
 175                     fprintf(out, "@ subtitle \"%s\"\n", subtitle);
 176                 }
 177             }
 178             if (g == ngraphs - 1)
 179             {
 180                 fprintf(out, "@ xaxis  label \"%s\"\n", xlabel.c_str());
 181             }
 182             else
 183             {
 184                 fprintf(out, "@ xaxis  ticklabel off\n");
 185             }
 186             if (n > 1)
 187             {
 188                 fprintf(out, "@ world xmin %g\n", x[0] * scale_x);
 189                 fprintf(out, "@ world xmax %g\n", x[n - 1] * scale_x);
 190                 fprintf(out, "@ world ymin %g\n", ymin);
 191                 fprintf(out, "@ world ymax %g\n", ymax);
 192             }
 193             fprintf(out, "@ view xmin 0.15\n");
 194             fprintf(out, "@ view xmax 0.85\n");
 195             fprintf(out, "@ view ymin %g\n", 0.15 + (ngraphs - 1 - g) * 0.7 / ngraphs);
 196             fprintf(out, "@ view ymax %g\n", 0.15 + (ngraphs - g) * 0.7 / ngraphs);
 197             fprintf(out, "@ yaxis  label \"%s\"\n", ylabel[g]);
 198             fprintf(out, "@ xaxis tick major %g\n", xsp);
 199             fprintf(out, "@ xaxis tick minor %g\n", xsp / 2);
 200             fprintf(out, "@ xaxis ticklabel start type spec\n");
 201             fprintf(out, "@ xaxis ticklabel start %g\n", std::ceil(ymin / xsp) * xsp);
 202             fprintf(out, "@ yaxis tick major %g\n", ysp);
 203             fprintf(out, "@ yaxis tick minor %g\n", ysp / 2);
 204             fprintf(out, "@ yaxis ticklabel start type spec\n");
 205             fprintf(out, "@ yaxis ticklabel start %g\n", std::ceil(ymin / ysp) * ysp);
 206             if ((ymin < 0) && (ymax > 0))
 207             {
 208                 fprintf(out, "@ zeroxaxis bar on\n");
 209                 fprintf(out, "@ zeroxaxis bar linestyle 3\n");
 210             }
 211         }
 212         for (s = 0; s < nsetspergraph; s++)
 213         {
 214             for (i = 0; i < n; i++)
 215             {
 216                 if (bSplit && i > 0 && std::abs(x[i]) < 1e-5)
 217                 {
 218                     fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 219                 }
 220                 fprintf(out, "%10.4f %10.5f\n", x[i] * scale_x, y ? y[g][i] : sy[g][s][i]);
 221             }
 222             fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 223         }
 224     }
 225     gmx_ffclose(out);
 226 }
 227
 228 static void
 229 compare(int natoms, int n1, rvec** eigvec1, int n2, rvec** eigvec2, real* eigval1, int neig1, real* eigval2, int neig2)
 230 {
 231     int    n;
 232     int    i, j, k;
 233     double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
 234
 235     n = std::min(n1, n2);
 236
 237     n = std::min(n, std::min(neig1, neig2));
 238     fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
 239
 240     sum1 = 0;
 241     for (i = 0; i < n; i++)
 242     {
 243         if (eigval1[i] < 0)
 244         {
 245             eigval1[i] = 0;
 246         }
 247         sum1 += eigval1[i];
 248         eigval1[i] = std::sqrt(eigval1[i]);
 249     }
 250     trace1 = sum1;
 251     for (i = n; i < neig1; i++)
 252     {
 253         trace1 += eigval1[i];
 254     }
 255     sum2 = 0;
 256     for (i = 0; i < n; i++)
 257     {
 258         if (eigval2[i] < 0)
 259         {
 260             eigval2[i] = 0;
 261         }
 262         sum2 += eigval2[i];
 263         eigval2[i] = std::sqrt(eigval2[i]);
 264     }
 265     trace2 = sum2;
 266
 267
 268     // The static analyzer appears to be confused by the fact that the loop below
 269     // starts from n instead of 0. However, given all the complex code it's
 270     // better to be safe than sorry, so we check it with an assert.
 271     // If we are in this comparison routine in the first place, neig2 should not be 0,
 272     // so eigval2 should always be a valid pointer.
 273     GMX_RELEASE_ASSERT(eigval2 != nullptr, "NULL pointer provided for eigval2");
 274
 275     for (i = n; i < neig2; i++)
 276     {
 277         trace2 += eigval2[i];
 278     }
 279
 280     fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
 281     if (neig1 != n || neig2 != n)
 282     {
 283         fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
 284                 gmx::roundToInt(100 * sum1 / trace1), gmx::roundToInt(100 * sum2 / trace2));
 285     }
 286     fprintf(stdout, "Square root of the traces: %g and %g\n", std::sqrt(sum1), std::sqrt(sum2));
 287
 288     sab = 0;
 289     for (i = 0; i < n; i++)
 290     {
 291         tmp = 0;
 292         for (j = 0; j < n; j++)
 293         {
 294             ip = 0;
 295             for (k = 0; k < natoms; k++)
 296             {
 297                 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
 298             }
 299             tmp += eigval2[j] * ip * ip;
 300         }
 301         sab += eigval1[i] * tmp;
 302     }
 303
 304     samsb2 = sum1 + sum2 - 2 * sab;
 305     if (samsb2 < 0)
 306     {
 307         samsb2 = 0;
 308     }
 309
 310     fprintf(stdout, "The overlap of the covariance matrices:\n");
 311     fprintf(stdout, "  normalized:  %.3f\n", 1 - std::sqrt(samsb2 / (sum1 + sum2)));
 312     tmp = 1 - sab / std::sqrt(sum1 * sum2);
 313     if (tmp < 0)
 314     {
 315         tmp = 0;
 316     }
 317     fprintf(stdout, "       shape:  %.3f\n", 1 - std::sqrt(tmp));
 318 }
 319
 320
 321 static void inprod_matrix(const char* matfile,
 322                           int         natoms,
 323                           int         nvec1,
 324                           int*        eignr1,
 325                           rvec**      eigvec1,
 326                           int         nvec2,
 327                           const int*  eignr2,
 328                           rvec**      eigvec2,
 329                           gmx_bool    bSelect,
 330                           int         noutvec,
 331                           const int*  outvec)
 332 {
 333     FILE*  out;
 334     real** mat;
 335     int    i, x1, y1, x, y, nlevels;
 336     int    nx, ny;
 337     real   inp, *t_x, *t_y, maxval;
 338     t_rgb  rlo, rhi;
 339
 340     snew(t_y, nvec2);
 341     if (bSelect)
 342     {
 343         nx = noutvec;
 344         ny = 0;
 345         for (y1 = 0; y1 < nx; y1++)
 346         {
 347             if (outvec[y1] < nvec2)
 348             {
 349                 t_y[ny] = eignr2[outvec[y1]] + 1;
 350                 ny++;
 351             }
 352         }
 353     }
 354     else
 355     {
 356         nx = nvec1;
 357         ny = nvec2;
 358         for (y = 0; y < ny; y++)
 359         {
 360             t_y[y] = eignr2[y] + 1;
 361         }
 362     }
 363
 364     fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n", nx, nvec2);
 365
 366     snew(mat, nx);
 367     snew(t_x, nx);
 368     maxval = 0;
 369     for (x1 = 0; x1 < nx; x1++)
 370     {
 371         snew(mat[x1], ny);
 372         if (bSelect)
 373         {
 374             x = outvec[x1];
 375         }
 376         else
 377         {
 378             x = x1;
 379         }
 380         t_x[x1] = eignr1[x] + 1;
 381         fprintf(stderr, " %d", eignr1[x] + 1);
 382         for (y1 = 0; y1 < ny; y1++)
 383         {
 384             inp = 0;
 385             if (bSelect)
 386             {
 387                 while (outvec[y1] >= nvec2)
 388                 {
 389                     y1++;
 390                 }
 391                 y = outvec[y1];
 392             }
 393             else
 394             {
 395                 y = y1;
 396             }
 397             for (i = 0; i < natoms; i++)
 398             {
 399                 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
 400             }
 401             mat[x1][y1] = std::abs(inp);
 402             if (mat[x1][y1] > maxval)
 403             {
 404                 maxval = mat[x1][y1];
 405             }
 406         }
 407     }
 408     fprintf(stderr, "\n");
 409     rlo.r   = 1;
 410     rlo.g   = 1;
 411     rlo.b   = 1;
 412     rhi.r   = 0;
 413     rhi.g   = 0;
 414     rhi.b   = 0;
 415     nlevels = 41;
 416     out     = gmx_ffopen(matfile, "w");
 417     write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2", nx, ny, t_x, t_y,
 418               mat, 0.0, maxval, rlo, rhi, &nlevels);
 419     gmx_ffclose(out);
 420 }
 421
 422 static void overlap(const char*             outfile,
 423                     int                     natoms,
 424                     rvec**                  eigvec1,
 425                     int                     nvec2,
 426                     int*                    eignr2,
 427                     rvec**                  eigvec2,
 428                     int                     noutvec,
 429                     int*                    outvec,
 430                     const gmx_output_env_t* oenv)
 431 {
 432     FILE* out;
 433     int   i, v, vec, x;
 434     real  overlap, inp;
 435
 436     fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
 437     for (i = 0; i < noutvec; i++)
 438     {
 439         fprintf(stderr, "%d ", outvec[i] + 1);
 440     }
 441     fprintf(stderr, "\n");
 442
 443     out = xvgropen(outfile, "Subspace overlap", "Eigenvectors of trajectory 2", "Overlap", oenv);
 444     if (output_env_get_print_xvgr_codes(oenv))
 445     {
 446         fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
 447     }
 448     overlap = 0;
 449     for (x = 0; x < nvec2; x++)
 450     {
 451         for (v = 0; v < noutvec; v++)
 452         {
 453             vec = outvec[v];
 454             inp = 0;
 455             for (i = 0; i < natoms; i++)
 456             {
 457                 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
 458             }
 459             overlap += gmx::square(inp);
 460         }
 461         fprintf(out, "%5d  %5.3f\n", eignr2[x] + 1, overlap / noutvec);
 462     }
 463
 464     xvgrclose(out);
 465 }
 466
 467 static void project(const char*             trajfile,
 468                     const t_topology*       top,
 469                     int                     ePBC,
 470                     matrix                  topbox,
 471                     const char*             projfile,
 472                     const char*             twodplotfile,
 473                     const char*             threedplotfile,
 474                     const char*             filterfile,
 475                     int                     skip,
 476                     const char*             extremefile,
 477                     gmx_bool                bExtrAll,
 478                     real                    extreme,
 479                     int                     nextr,
 480                     const t_atoms*          atoms,
 481                     int                     natoms,
 482                     int*                    index,
 483                     gmx_bool                bFit,
 484                     rvec*                   xref,
 485                     int                     nfit,
 486                     int*                    ifit,
 487                     real*                   w_rls,
 488                     const real*             sqrtm,
 489                     rvec*                   xav,
 490                     int*                    eignr,
 491                     rvec**                  eigvec,
 492                     int                     noutvec,
 493                     int*                    outvec,
 494                     gmx_bool                bSplit,
 495                     const gmx_output_env_t* oenv)
 496 {
 497     FILE*        xvgrout = nullptr;
 498     int          nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
 499     t_trxstatus* out = nullptr;
 500     t_trxstatus* status;
 501     int          noutvec_extr, imin, imax;
 502     real *       pmin, *pmax;
 503     int*         all_at;
 504     matrix       box;
 505     rvec *       xread, *x;
 506     real         t, inp, **inprod = nullptr;
 507     char         str[STRLEN], str2[STRLEN], *c;
 508     const char** ylabel;
 509     real         fact;
 510     gmx_rmpbc_t  gpbc = nullptr;
 511
 512     snew(x, natoms);
 513
 514     if (bExtrAll)
 515     {
 516         noutvec_extr = noutvec;
 517     }
 518     else
 519     {
 520         noutvec_extr = 1;
 521     }
 522
 523
 524     if (trajfile)
 525     {
 526         snew(inprod, noutvec + 1);
 527
 528         if (filterfile)
 529         {
 530             fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n", filterfile);
 531             for (i = 0; i < noutvec; i++)
 532             {
 533                 fprintf(stderr, "%d ", outvec[i] + 1);
 534             }
 535             fprintf(stderr, "\n");
 536             out = open_trx(filterfile, "w");
 537         }
 538         snew_size = 0;
 539         nfr       = 0;
 540         nframes   = 0;
 541         nat       = read_first_x(oenv, &status, trajfile, &t, &xread, box);
 542         if (nat > atoms->nr)
 543         {
 544             gmx_fatal(FARGS,
 545                       "the number of atoms in your trajectory (%d) is larger than the number of "
 546                       "atoms in your structure file (%d)",
 547                       nat, atoms->nr);
 548         }
 549         snew(all_at, nat);
 550
 551         if (top)
 552         {
 553             gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
 554         }
 555
 556         for (i = 0; i < nat; i++)
 557         {
 558             all_at[i] = i;
 559         }
 560         do
 561         {
 562             if (nfr % skip == 0)
 563             {
 564                 if (top)
 565                 {
 566                     gmx_rmpbc(gpbc, nat, box, xread);
 567                 }
 568                 if (nframes >= snew_size)
 569                 {
 570                     snew_size += 100;
 571                     for (i = 0; i < noutvec + 1; i++)
 572                     {
 573                         srenew(inprod[i], snew_size);
 574                     }
 575                 }
 576                 inprod[noutvec][nframes] = t;
 577                 /* calculate x: a fitted struture of the selected atoms */
 578                 if (bFit)
 579                 {
 580                     reset_x(nfit, ifit, nat, nullptr, xread, w_rls);
 581                     do_fit(nat, w_rls, xref, xread);
 582                 }
 583                 for (i = 0; i < natoms; i++)
 584                 {
 585                     copy_rvec(xread[index[i]], x[i]);
 586                 }
 587
 588                 for (v = 0; v < noutvec; v++)
 589                 {
 590                     vec = outvec[v];
 591                     /* calculate (mass-weighted) projection */
 592                     inp = 0;
 593                     for (i = 0; i < natoms; i++)
 594                     {
 595                         inp += (eigvec[vec][i][0] * (x[i][0] - xav[i][0])
 596                                 + eigvec[vec][i][1] * (x[i][1] - xav[i][1])
 597                                 + eigvec[vec][i][2] * (x[i][2] - xav[i][2]))
 598                                * sqrtm[i];
 599                     }
 600                     inprod[v][nframes] = inp;
 601                 }
 602                 if (filterfile)
 603                 {
 604                     for (i = 0; i < natoms; i++)
 605                     {
 606                         for (d = 0; d < DIM; d++)
 607                         {
 608                             /* misuse xread for output */
 609                             xread[index[i]][d] = xav[i][d];
 610                             for (v = 0; v < noutvec; v++)
 611                             {
 612                                 xread[index[i]][d] +=
 613                                         inprod[v][nframes] * eigvec[outvec[v]][i][d] / sqrtm[i];
 614                             }
 615                         }
 616                     }
 617                     write_trx(out, natoms, index, atoms, 0, t, box, xread, nullptr, nullptr);
 618                 }
 619                 nframes++;
 620             }
 621             nfr++;
 622         } while (read_next_x(oenv, status, &t, xread, box));
 623         close_trx(status);
 624         sfree(x);
 625         if (filterfile)
 626         {
 627             close_trx(out);
 628         }
 629     }
 630     else
 631     {
 632         snew(xread, atoms->nr);
 633     }
 634
 635     if (top)
 636     {
 637         gmx_rmpbc_done(gpbc);
 638     }
 639
 640
 641     if (projfile)
 642     {
 643         GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL if projfile is non-NULL");
 644         snew(ylabel, noutvec);
 645         for (v = 0; v < noutvec; v++)
 646         {
 647             sprintf(str, "vec %d", eignr[outvec[v]] + 1);
 648             ylabel[v] = gmx_strdup(str);
 649         }
 650         sprintf(str, "projection on eigenvectors (%s)", proj_unit);
 651         write_xvgr_graphs(projfile, noutvec, 1, str, nullptr, output_env_get_xvgr_tlabel(oenv),
 652                           ylabel, nframes, inprod[noutvec], inprod, nullptr,
 653                           output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
 654     }
 655
 656     if (twodplotfile)
 657     {
 658         sprintf(str, "projection on eigenvector %d (%s)", eignr[outvec[0]] + 1, proj_unit);
 659         sprintf(str2, "projection on eigenvector %d (%s)", eignr[outvec[noutvec - 1]] + 1, proj_unit);
 660         xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2, oenv);
 661         for (i = 0; i < nframes; i++)
 662         {
 663             if (bSplit && i > 0 && std::abs(inprod[noutvec][i]) < 1e-5)
 664             {
 665                 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 666             }
 667             fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec - 1][i]);
 668         }
 669         xvgrclose(xvgrout);
 670     }
 671
 672     if (threedplotfile)
 673     {
 674         t_atoms  atoms;
 675         rvec*    x;
 676         real*    b = nullptr;
 677         matrix   box;
 678         char *   resnm, *atnm;
 679         gmx_bool bPDB, b4D;
 680         FILE*    out;
 681
 682         if (noutvec < 3)
 683         {
 684             gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
 685         }
 686
 687         /* initialize */
 688         bPDB = fn2ftp(threedplotfile) == efPDB;
 689         clear_mat(box);
 690         box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
 691
 692         b4D = bPDB && (noutvec >= 4);
 693         if (b4D)
 694         {
 695             fprintf(stderr,
 696                     "You have selected four or more eigenvectors:\n"
 697                     "fourth eigenvector will be plotted "
 698                     "in bfactor field of pdb file\n");
 699             sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d", eignr[outvec[0]] + 1,
 700                     eignr[outvec[1]] + 1, eignr[outvec[2]] + 1, eignr[outvec[3]] + 1);
 701         }
 702         else
 703         {
 704             sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d", eignr[outvec[0]] + 1,
 705                     eignr[outvec[1]] + 1, eignr[outvec[2]] + 1);
 706         }
 707         init_t_atoms(&atoms, nframes, FALSE);
 708         snew(x, nframes);
 709         snew(b, nframes);
 710         atnm  = gmx_strdup("C");
 711         resnm = gmx_strdup("PRJ");
 712
 713         if (nframes > 10000)
 714         {
 715             fact = 10000.0 / nframes;
 716         }
 717         else
 718         {
 719             fact = 1.0;
 720         }
 721
 722         for (i = 0; i < nframes; i++)
 723         {
 724             atoms.atomname[i]     = &atnm;
 725             atoms.atom[i].resind  = i;
 726             atoms.resinfo[i].name = &resnm;
 727             atoms.resinfo[i].nr   = static_cast<int>(std::ceil(i * fact));
 728             atoms.resinfo[i].ic   = ' ';
 729             x[i][XX]              = inprod[0][i];
 730             x[i][YY]              = inprod[1][i];
 731             x[i][ZZ]              = inprod[2][i];
 732             if (b4D)
 733             {
 734                 b[i] = inprod[3][i];
 735             }
 736         }
 737         if ((b4D || bSplit) && bPDB)
 738         {
 739             GMX_RELEASE_ASSERT(inprod != nullptr,
 740                                "inprod must be non-NULL with 4D or split PDB output options");
 741
 742             out = gmx_ffopen(threedplotfile, "w");
 743             fprintf(out, "HEADER    %s\n", str);
 744             if (b4D)
 745             {
 746                 fprintf(out, "REMARK    %s\n", "fourth dimension plotted as B-factor");
 747             }
 748             j = 0;
 749             for (i = 0; i < atoms.nr; i++)
 750             {
 751                 if (j > 0 && bSplit && std::abs(inprod[noutvec][i]) < 1e-5)
 752                 {
 753                     fprintf(out, "TER\n");
 754                     j = 0;
 755                 }
 756                 gmx_fprintf_pdb_atomline(out, epdbATOM, i + 1, "C", ' ', "PRJ", ' ', j + 1, ' ',
 757                                          10 * x[i][XX], 10 * x[i][YY], 10 * x[i][ZZ], 1.0,
 758                                          10 * b[i], "");
 759                 if (j > 0)
 760                 {
 761                     fprintf(out, "CONECT%5d%5d\n", i, i + 1);
 762                 }
 763                 j++;
 764             }
 765             fprintf(out, "TER\n");
 766             gmx_ffclose(out);
 767         }
 768         else
 769         {
 770             write_sto_conf(threedplotfile, str, &atoms, x, nullptr, ePBC, box);
 771         }
 772         done_atom(&atoms);
 773     }
 774
 775     if (extremefile)
 776     {
 777         snew(pmin, noutvec_extr);
 778         snew(pmax, noutvec_extr);
 779         if (extreme == 0)
 780         {
 781             GMX_RELEASE_ASSERT(inprod != nullptr, "inprod must be non-NULL");
 782             fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
 783             fprintf(stderr, "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
 784             imin = 0;
 785             imax = 0;
 786             for (v = 0; v < noutvec_extr; v++)
 787             {
 788                 for (i = 0; i < nframes; i++)
 789                 {
 790                     if (inprod[v][i] < inprod[v][imin])
 791                     {
 792                         imin = i;
 793                     }
 794                     if (inprod[v][i] > inprod[v][imax])
 795                     {
 796                         imax = i;
 797                     }
 798                 }
 799                 pmin[v] = inprod[v][imin];
 800                 pmax[v] = inprod[v][imax];
 801                 fprintf(stderr, "%7d     %10.6f %10d %10.6f %10d\n", eignr[outvec[v]] + 1, pmin[v],
 802                         imin, pmax[v], imax);
 803             }
 804         }
 805         else
 806         {
 807             pmin[0] = -extreme;
 808             pmax[0] = extreme;
 809         }
 810         /* build format string for filename: */
 811         std::strcpy(str, extremefile); /* copy filename */
 812         c = std::strrchr(str, '.');    /* find where extention begins */
 813         std::strcpy(str2, c);          /* get extention */
 814         sprintf(c, "%%d%s", str2);     /* append '%s' and extention to filename */
 815         for (v = 0; v < noutvec_extr; v++)
 816         {
 817             /* make filename using format string */
 818             if (noutvec_extr == 1)
 819             {
 820                 std::strcpy(str2, extremefile);
 821             }
 822             else
 823             {
 824                 sprintf(str2, str, eignr[outvec[v]] + 1);
 825             }
 826             fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n", nextr, outvec[v] + 1, str2);
 827             out = open_trx(str2, "w");
 828             for (frame = 0; frame < nextr; frame++)
 829             {
 830                 if ((extreme == 0) && (nextr <= 3))
 831                 {
 832                     for (i = 0; i < natoms; i++)
 833                     {
 834                         atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
 835                     }
 836                 }
 837                 for (i = 0; i < natoms; i++)
 838                 {
 839                     for (d = 0; d < DIM; d++)
 840                     {
 841                         xread[index[i]][d] =
 842                                 (xav[i][d]
 843                                  + (pmin[v] * (nextr - frame - 1) + pmax[v] * frame) / (nextr - 1)
 844                                            * eigvec[outvec[v]][i][d] / sqrtm[i]);
 845                     }
 846                 }
 847                 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, nullptr, nullptr);
 848             }
 849             close_trx(out);
 850         }
 851         sfree(pmin);
 852         sfree(pmax);
 853     }
 854     fprintf(stderr, "\n");
 855 }
 856
 857 static void components(const char*             outfile,
 858                        int                     natoms,
 859                        int*                    eignr,
 860                        rvec**                  eigvec,
 861                        int                     noutvec,
 862                        const int*              outvec,
 863                        const gmx_output_env_t* oenv)
 864 {
 865     int          g, s, v, i;
 866     real *       x, ***y;
 867     char         str[STRLEN];
 868     const char** ylabel;
 869
 870     fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
 871
 872     snew(ylabel, noutvec);
 873     snew(y, noutvec);
 874     snew(x, natoms);
 875     for (i = 0; i < natoms; i++)
 876     {
 877         x[i] = i + 1;
 878     }
 879     for (g = 0; g < noutvec; g++)
 880     {
 881         v = outvec[g];
 882         sprintf(str, "vec %d", eignr[v] + 1);
 883         ylabel[g] = gmx_strdup(str);
 884         snew(y[g], 4);
 885         for (s = 0; s < 4; s++)
 886         {
 887             snew(y[g][s], natoms);
 888         }
 889         for (i = 0; i < natoms; i++)
 890         {
 891             y[g][0][i] = norm(eigvec[v][i]);
 892             for (s = 0; s < 3; s++)
 893             {
 894                 y[g][s + 1][i] = eigvec[v][i][s];
 895             }
 896         }
 897     }
 898     write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
 899                       "black: total, red: x, green: y, blue: z", "Atom number", ylabel, natoms, x,
 900                       nullptr, y, 1, FALSE, FALSE, oenv);
 901     fprintf(stderr, "\n");
 902 }
 903
 904 static void rmsf(const char*             outfile,
 905                  int                     natoms,
 906                  const real*             sqrtm,
 907                  int*                    eignr,
 908                  rvec**                  eigvec,
 909                  int                     noutvec,
 910                  const int*              outvec,
 911                  real*                   eigval,
 912                  int                     neig,
 913                  const gmx_output_env_t* oenv)
 914 {
 915     int          g, v, i;
 916     real *       x, **y;
 917     char         str[STRLEN];
 918     const char** ylabel;
 919
 920     for (i = 0; i < neig; i++)
 921     {
 922         if (eigval[i] < 0)
 923         {
 924             eigval[i] = 0;
 925         }
 926     }
 927
 928     fprintf(stderr, "Writing rmsf to %s\n", outfile);
 929
 930     snew(ylabel, noutvec);
 931     snew(y, noutvec);
 932     snew(x, natoms);
 933     for (i = 0; i < natoms; i++)
 934     {
 935         x[i] = i + 1;
 936     }
 937     for (g = 0; g < noutvec; g++)
 938     {
 939         v = outvec[g];
 940         if (eignr[v] >= neig)
 941         {
 942             gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)",
 943                       eignr[v] + 1, neig);
 944         }
 945         sprintf(str, "vec %d", eignr[v] + 1);
 946         ylabel[g] = gmx_strdup(str);
 947         snew(y[g], natoms);
 948         for (i = 0; i < natoms; i++)
 949         {
 950             y[g][i] = std::sqrt(eigval[eignr[v]] * norm2(eigvec[v][i])) / sqrtm[i];
 951         }
 952     }
 953     write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", nullptr, "Atom number", ylabel,
 954                       natoms, x, y, nullptr, 1, TRUE, FALSE, oenv);
 955     fprintf(stderr, "\n");
 956 }
 957
 958 int gmx_anaeig(int argc, char* argv[])
 959 {
 960     static const char* desc[] = {
 961         "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
 962         "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
 963         "([gmx-nmeig]).[PAR]",
 964
 965         "When a trajectory is projected on eigenvectors, all structures are",
 966         "fitted to the structure in the eigenvector file, if present, otherwise",
 967         "to the structure in the structure file. When no run input file is",
 968         "supplied, periodicity will not be taken into account. Most analyses",
 969         "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
 970         "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
 971
 972         "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
 973         "[TT]-first[tt] to [TT]-last[tt].[PAR]",
 974
 975         "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
 976         "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
 977
 978         "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
 979         "[TT]-first[tt] to [TT]-last[tt].",
 980         "The projections of a trajectory on the eigenvectors of its",
 981         "covariance matrix are called principal components (pc's).",
 982         "It is often useful to check the cosine content of the pc's,",
 983         "since the pc's of random diffusion are cosines with the number",
 984         "of periods equal to half the pc index.",
 985         "The cosine content of the pc's can be calculated with the program",
 986         "[gmx-analyze].[PAR]",
 987
 988         "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
 989         "[TT]-first[tt] and [TT]-last[tt].[PAR]",
 990
 991         "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
 992         "three selected eigenvectors.[PAR]",
 993
 994         "[TT]-filt[tt]: filter the trajectory to show only the motion along",
 995         "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
 996
 997         "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
 998         "on the average structure and interpolate [TT]-nframes[tt] frames",
 999         "between them, or set your own extremes with [TT]-max[tt]. The",
1000         "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
1001         "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1002         "will be written to separate files. Chain identifiers will be added",
1003         "when writing a [REF].pdb[ref] file with two or three structures (you",
1004         "can use [TT]rasmol -nmrpdb[tt] to view such a [REF].pdb[ref] file).[PAR]",
1005
1006         "Overlap calculations between covariance analysis",
1007         "^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^",
1008         "",
1009         "[BB]Note:[bb] the analysis should use the same fitting structure",
1010         "",
1011         "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1012         "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1013         "in file [TT]-v[tt].[PAR]",
1014
1015         "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1016         "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1017         "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1018         "have been set explicitly.[PAR]",
1019
1020         "When [TT]-v[tt] and [TT]-v2[tt] are given, a single number for the",
1021         "overlap between the covariance matrices is generated. Note that the",
1022         "eigenvalues are by default read from the timestamp field in the",
1023         "eigenvector input files, but when [TT]-eig[tt], or [TT]-eig2[tt] are",
1024         "given, the corresponding eigenvalues are used instead. The formulas are::",
1025         "",
1026         "         difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))",
1027         " normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))",
1028         "      shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))",
1029         "",
1030         "where M1 and M2 are the two covariance matrices and tr is the trace",
1031         "of a matrix. The numbers are proportional to the overlap of the square",
1032         "root of the fluctuations. The normalized overlap is the most useful",
1033         "number, it is 1 for identical matrices and 0 when the sampled",
1034         "subspaces are orthogonal.[PAR]",
1035         "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1036         "computed based on the Quasiharmonic approach and based on",
1037         "Schlitter's formula."
1038     };
1039     static int      first = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1040     static real     max = 0.0, temp = 298.15;
1041     static gmx_bool bSplit = FALSE, bEntropy = FALSE;
1042     t_pargs         pa[] = {
1043         { "-first", FALSE, etINT, { &first }, "First eigenvector for analysis (-1 is select)" },
1044         { "-last", FALSE, etINT, { &last }, "Last eigenvector for analysis (-1 is till the last)" },
1045         { "-skip", FALSE, etINT, { &skip }, "Only analyse every nr-th frame" },
1046         { "-max",
1047           FALSE,
1048           etREAL,
1049           { &max },
1050           "Maximum for projection of the eigenvector on the average structure, "
1051           "max=0 gives the extremes" },
1052         { "-nframes", FALSE, etINT, { &nextr }, "Number of frames for the extremes output" },
1053         { "-split", FALSE, etBOOL, { &bSplit }, "Split eigenvector projections where time is zero" },
1054         { "-entropy",
1055           FALSE,
1056           etBOOL,
1057           { &bEntropy },
1058           "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1059         { "-temp", FALSE, etREAL, { &temp }, "Temperature for entropy calculations" },
1060         { "-nevskip",
1061           FALSE,
1062           etINT,
1063           { &nskip },
1064           "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic "
1065           "approximation. When you do a rotational and/or translational fit prior to the "
1066           "covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which "
1067           "should not be taken into account when computing the entropy." }
1068     };
1069 #define NPA asize(pa)
1070
1071     t_topology        top;
1072     int               ePBC  = -1;
1073     const t_atoms*    atoms = nullptr;
1074     rvec *            xtop, *xref1, *xref2, *xrefp = nullptr;
1075     gmx_bool          bDMR1, bDMA1, bDMR2, bDMA2;
1076     int               nvec1, nvec2, *eignr1 = nullptr, *eignr2 = nullptr;
1077     rvec *            xav1, *xav2, **eigvec1 = nullptr, **eigvec2 = nullptr;
1078     matrix            topbox;
1079     real              totmass, *sqrtm, *w_rls, t;
1080     int               natoms;
1081     char*             grpname;
1082     const char*       indexfile;
1083     int               i, j, d;
1084     int               nout, *iout, noutvec, *outvec, nfit;
1085     int *             index = nullptr, *ifit = nullptr;
1086     const char *      VecFile, *Vec2File, *topfile;
1087     const char *      EigFile, *Eig2File;
1088     const char *      CompFile, *RmsfFile, *ProjOnVecFile;
1089     const char *      TwoDPlotFile, *ThreeDPlotFile;
1090     const char *      FilterFile, *ExtremeFile;
1091     const char *      OverlapFile, *InpMatFile;
1092     gmx_bool          bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1093     gmx_bool          bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1094     real *            eigval1 = nullptr, *eigval2 = nullptr;
1095     int               neig1, neig2;
1096     double**          xvgdata;
1097     gmx_output_env_t* oenv;
1098     gmx_rmpbc_t       gpbc;
1099
1100     t_filenm fnm[] = {
1101         { efTRN, "-v", "eigenvec", ffREAD },      { efTRN, "-v2", "eigenvec2", ffOPTRD },
1102         { efTRX, "-f", nullptr, ffOPTRD },        { efTPS, nullptr, nullptr, ffOPTRD },
1103         { efNDX, nullptr, nullptr, ffOPTRD },     { efXVG, "-eig", "eigenval", ffOPTRD },
1104         { efXVG, "-eig2", "eigenval2", ffOPTRD }, { efXVG, "-comp", "eigcomp", ffOPTWR },
1105         { efXVG, "-rmsf", "eigrmsf", ffOPTWR },   { efXVG, "-proj", "proj", ffOPTWR },
1106         { efXVG, "-2d", "2dproj", ffOPTWR },      { efSTO, "-3d", "3dproj.pdb", ffOPTWR },
1107         { efTRX, "-filt", "filtered", ffOPTWR },  { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1108         { efXVG, "-over", "overlap", ffOPTWR },   { efXPM, "-inpr", "inprod", ffOPTWR }
1109     };
1110 #define NFILE asize(fnm)
1111
1112     if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW, NFILE, fnm,
1113                            NPA, pa, asize(desc), desc, 0, nullptr, &oenv))
1114     {
1115         return 0;
1116     }
1117
1118     indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1119
1120     VecFile        = opt2fn("-v", NFILE, fnm);
1121     Vec2File       = opt2fn_null("-v2", NFILE, fnm);
1122     topfile        = ftp2fn(efTPS, NFILE, fnm);
1123     EigFile        = opt2fn_null("-eig", NFILE, fnm);
1124     Eig2File       = opt2fn_null("-eig2", NFILE, fnm);
1125     CompFile       = opt2fn_null("-comp", NFILE, fnm);
1126     RmsfFile       = opt2fn_null("-rmsf", NFILE, fnm);
1127     ProjOnVecFile  = opt2fn_null("-proj", NFILE, fnm);
1128     TwoDPlotFile   = opt2fn_null("-2d", NFILE, fnm);
1129     ThreeDPlotFile = opt2fn_null("-3d", NFILE, fnm);
1130     FilterFile     = opt2fn_null("-filt", NFILE, fnm);
1131     ExtremeFile    = opt2fn_null("-extr", NFILE, fnm);
1132     OverlapFile    = opt2fn_null("-over", NFILE, fnm);
1133     InpMatFile     = ftp2fn_null(efXPM, NFILE, fnm);
1134
1135     bProj = (ProjOnVecFile != nullptr) || (TwoDPlotFile != nullptr) || (ThreeDPlotFile != nullptr)
1136             || (FilterFile != nullptr) || (ExtremeFile != nullptr);
1137     bFirstLastSet = opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1138     bFirstToLast  = (CompFile != nullptr) || (RmsfFile != nullptr) || (ProjOnVecFile != nullptr)
1139                    || (FilterFile != nullptr) || (OverlapFile != nullptr)
1140                    || (((ExtremeFile != nullptr) || (InpMatFile != nullptr)) && bFirstLastSet);
1141     bVec2 = (Vec2File != nullptr) || (OverlapFile != nullptr) || (InpMatFile != nullptr);
1142     bM    = (RmsfFile != nullptr) || bProj;
1143     bTraj = (ProjOnVecFile != nullptr) || (FilterFile != nullptr)
1144             || ((ExtremeFile != nullptr) && (max == 0)) || (TwoDPlotFile != nullptr)
1145             || (ThreeDPlotFile != nullptr);
1146     bIndex = bM || bProj;
1147     bTPS   = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj || (FilterFile != nullptr)
1148            || (bIndex && (indexfile != nullptr));
1149     bCompare = (Vec2File != nullptr) || (Eig2File != nullptr);
1150     bPDB3D   = fn2ftp(ThreeDPlotFile) == efPDB;
1151
1152     read_eigenvectors(VecFile, &natoms, &bFit1, &xref1, &bDMR1, &xav1, &bDMA1, &nvec1, &eignr1,
1153                       &eigvec1, &eigval1);
1154     neig1 = std::min(nvec1, DIM * natoms);
1155     if (nvec1 != DIM * natoms)
1156     {
1157         fprintf(stderr,
1158                 "Warning: number of eigenvectors %d does not match three times\n"
1159                 "the number of atoms %d in %s. Using %d eigenvectors.\n\n",
1160                 nvec1, natoms, VecFile, neig1);
1161     }
1162
1163     /* Overwrite eigenvalues from separate files if the user provides them */
1164     if (EigFile != nullptr)
1165     {
1166         int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1167         if (neig_tmp != neig1)
1168         {
1169             fprintf(stderr,
1170                     "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n",
1171                     neig1, natoms);
1172         }
1173         neig1 = neig_tmp;
1174         srenew(eigval1, neig1);
1175         for (j = 0; j < neig1; j++)
1176         {
1177             real tmp   = eigval1[j];
1178             eigval1[j] = xvgdata[1][j];
1179             if (debug && (eigval1[j] != tmp))
1180             {
1181                 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n", j, tmp,
1182                         eigval1[j]);
1183             }
1184         }
1185         for (j = 0; j < i; j++)
1186         {
1187             sfree(xvgdata[j]);
1188         }
1189         sfree(xvgdata);
1190         fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1191     }
1192
1193     if (bEntropy)
1194     {
1195         if (bDMA1)
1196         {
1197             gmx_fatal(FARGS,
1198                       "Can not calculate entropies from mass-weighted eigenvalues, redo the "
1199                       "analysis without mass-weighting");
1200         }
1201         printf("The Entropy due to the Schlitter formula is %g J/mol K\n",
1202                calcSchlitterEntropy(gmx::arrayRefFromArray(eigval1, neig1), temp, FALSE));
1203         printf("The Entropy due to the Quasiharmonic analysis is %g J/mol K\n",
1204                calcQuasiHarmonicEntropy(gmx::arrayRefFromArray(eigval1, neig1), temp, FALSE, 1.0));
1205     }
1206
1207     if (bVec2)
1208     {
1209         if (!Vec2File)
1210         {
1211             gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1212         }
1213         int natoms2;
1214         read_eigenvectors(Vec2File, &natoms2, &bFit2, &xref2, &bDMR2, &xav2, &bDMA2, &nvec2,
1215                           &eignr2, &eigvec2, &eigval2);
1216
1217         neig2 = std::min(nvec2, DIM * natoms2);
1218         if (neig2 != neig1)
1219         {
1220             gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1221         }
1222     }
1223     else
1224     {
1225         nvec2 = 0;
1226         neig2 = 0;
1227     }
1228
1229     if (Eig2File != nullptr)
1230     {
1231         neig2 = read_xvg(Eig2File, &xvgdata, &i);
1232         srenew(eigval2, neig2);
1233         for (j = 0; j < neig2; j++)
1234         {
1235             eigval2[j] = xvgdata[1][j];
1236         }
1237         for (j = 0; j < i; j++)
1238         {
1239             sfree(xvgdata[j]);
1240         }
1241         sfree(xvgdata);
1242         fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1243     }
1244
1245
1246     if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1247     {
1248         bM = FALSE;
1249     }
1250     if ((xref1 == nullptr) && (bM || bTraj))
1251     {
1252         bTPS = TRUE;
1253     }
1254
1255     xtop  = nullptr;
1256     nfit  = 0;
1257     ifit  = nullptr;
1258     w_rls = nullptr;
1259
1260     if (!bTPS)
1261     {
1262         bTop = FALSE;
1263     }
1264     else
1265     {
1266         bTop  = read_tps_conf(ftp2fn(efTPS, NFILE, fnm), &top, &ePBC, &xtop, nullptr, topbox, bM);
1267         atoms = &top.atoms;
1268         gpbc  = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1269         gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1270         /* Fitting is only required for the projection */
1271         if (bProj && bFit1)
1272         {
1273             if (xref1 == nullptr)
1274             {
1275                 printf("\nNote: the structure in %s should be the same\n"
1276                        "      as the one used for the fit in g_covar\n",
1277                        topfile);
1278             }
1279             printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1280             get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1281
1282             snew(w_rls, atoms->nr);
1283             for (i = 0; (i < nfit); i++)
1284             {
1285                 if (bDMR1)
1286                 {
1287                     w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1288                 }
1289                 else
1290                 {
1291                     w_rls[ifit[i]] = 1.0;
1292                 }
1293             }
1294
1295             snew(xrefp, atoms->nr);
1296             if (xref1 != nullptr)
1297             {
1298                 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1299                 if (natoms != nfit)
1300                 {
1301                     gmx_fatal(FARGS,
1302                               "you selected a group with %d elements instead of %d, your selection "
1303                               "does not fit the reference structure in the eigenvector file.",
1304                               nfit, natoms);
1305                 }
1306                 for (i = 0; (i < nfit); i++)
1307                 {
1308                     copy_rvec(xref1[i], xrefp[ifit[i]]);
1309                 }
1310             }
1311             else
1312             {
1313                 /* The top coordinates are the fitting reference */
1314                 for (i = 0; (i < nfit); i++)
1315                 {
1316                     copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1317                 }
1318                 reset_x(nfit, ifit, atoms->nr, nullptr, xrefp, w_rls);
1319             }
1320         }
1321         gmx_rmpbc_done(gpbc);
1322     }
1323
1324     if (bIndex)
1325     {
1326         printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1327         get_index(atoms, indexfile, 1, &i, &index, &grpname);
1328         if (i != natoms)
1329         {
1330             gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1331         }
1332         printf("\n");
1333     }
1334
1335     snew(sqrtm, natoms);
1336     if (bM && bDMA1)
1337     {
1338         proj_unit = "u\\S1/2\\Nnm";
1339         for (i = 0; (i < natoms); i++)
1340         {
1341             sqrtm[i] = std::sqrt(atoms->atom[index[i]].m);
1342         }
1343     }
1344     else
1345     {
1346         proj_unit = "nm";
1347         for (i = 0; (i < natoms); i++)
1348         {
1349             sqrtm[i] = 1.0;
1350         }
1351     }
1352
1353     if (bVec2)
1354     {
1355         t       = 0;
1356         totmass = 0;
1357         for (i = 0; (i < natoms); i++)
1358         {
1359             for (d = 0; (d < DIM); d++)
1360             {
1361                 t += gmx::square((xav1[i][d] - xav2[i][d]) * sqrtm[i]);
1362                 totmass += gmx::square(sqrtm[i]);
1363             }
1364         }
1365         fprintf(stdout,
1366                 "RMSD (without fit) between the two average structures:"
1367                 " %.3f (nm)\n\n",
1368                 std::sqrt(t / totmass));
1369     }
1370
1371     if (last == -1)
1372     {
1373         last = natoms * DIM;
1374     }
1375     if (first > -1)
1376     {
1377         if (bFirstToLast)
1378         {
1379             /* make an index from first to last */
1380             nout = last - first + 1;
1381             snew(iout, nout);
1382             for (i = 0; i < nout; i++)
1383             {
1384                 iout[i] = first - 1 + i;
1385             }
1386         }
1387         else if (ThreeDPlotFile)
1388         {
1389             /* make an index of first+(0,1,2) and last */
1390             nout = bPDB3D ? 4 : 3;
1391             nout = std::min(last - first + 1, nout);
1392             snew(iout, nout);
1393             iout[0] = first - 1;
1394             iout[1] = first;
1395             if (nout > 3)
1396             {
1397                 iout[2] = first + 1;
1398             }
1399             iout[nout - 1] = last - 1;
1400         }
1401         else
1402         {
1403             /* make an index of first and last */
1404             nout = 2;
1405             snew(iout, nout);
1406             iout[0] = first - 1;
1407             iout[1] = last - 1;
1408         }
1409     }
1410     else
1411     {
1412         printf("Select eigenvectors for output, end your selection with 0\n");
1413         nout = -1;
1414         iout = nullptr;
1415
1416         do
1417         {
1418             nout++;
1419             srenew(iout, nout + 1);
1420             if (1 != scanf("%d", &iout[nout]))
1421             {
1422                 gmx_fatal(FARGS, "Error reading user input");
1423             }
1424             iout[nout]--;
1425         } while (iout[nout] >= 0);
1426
1427         printf("\n");
1428     }
1429     /* make an index of the eigenvectors which are present */
1430     snew(outvec, nout);
1431     noutvec = 0;
1432     for (i = 0; i < nout; i++)
1433     {
1434         j = 0;
1435         while ((j < nvec1) && (eignr1[j] != iout[i]))
1436         {
1437             j++;
1438         }
1439         if ((j < nvec1) && (eignr1[j] == iout[i]))
1440         {
1441             outvec[noutvec] = j;
1442             noutvec++;
1443         }
1444     }
1445     fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1446     if (noutvec <= 100)
1447     {
1448         fprintf(stderr, ":");
1449         for (j = 0; j < noutvec; j++)
1450         {
1451             fprintf(stderr, " %d", eignr1[outvec[j]] + 1);
1452         }
1453     }
1454     fprintf(stderr, "\n");
1455
1456     if (CompFile)
1457     {
1458         components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1459     }
1460
1461     if (RmsfFile)
1462     {
1463         rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1, neig1, oenv);
1464     }
1465
1466     if (bProj)
1467     {
1468         project(bTraj ? opt2fn("-f", NFILE, fnm) : nullptr, bTop ? &top : nullptr, ePBC, topbox,
1469                 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip, ExtremeFile,
1470                 bFirstLastSet, max, nextr, atoms, natoms, index, bFit1, xrefp, nfit, ifit, w_rls,
1471                 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit, oenv);
1472     }
1473
1474     if (OverlapFile)
1475     {
1476         overlap(OverlapFile, natoms, eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1477     }
1478
1479     if (InpMatFile)
1480     {
1481         inprod_matrix(InpMatFile, natoms, nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1482                       bFirstLastSet, noutvec, outvec);
1483     }
1484
1485     if (bCompare)
1486     {
1487         compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1488     }
1489
1490
1491     if (!CompFile && !bProj && !OverlapFile && !InpMatFile && !bCompare && !bEntropy)
1492     {
1493         fprintf(stderr,
1494                 "\nIf you want some output,"
1495                 " set one (or two or ...) of the output file options\n");
1496     }
1497
1498
1499     view_all(oenv, NFILE, fnm);
1500
1501     return 0;
1502 }