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