src/gromacs/gmxana/gmx_anaeig.c

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