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