src/gromacs/gmxana/gmx_anaeig.c

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