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