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