src/gromacs/gmxana/gmx_anaeig.cpp

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