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     if (!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         return 0;
1114     }
1115
1116     indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1117
1118     VecFile         = opt2fn("-v", NFILE, fnm);
1119     Vec2File        = opt2fn_null("-v2", NFILE, fnm);
1120     topfile         = ftp2fn(efTPS, NFILE, fnm);
1121     EigFile         = opt2fn_null("-eig", NFILE, fnm);
1122     Eig2File        = opt2fn_null("-eig2", NFILE, fnm);
1123     CompFile        = opt2fn_null("-comp", NFILE, fnm);
1124     RmsfFile        = opt2fn_null("-rmsf", NFILE, fnm);
1125     ProjOnVecFile   = opt2fn_null("-proj", NFILE, fnm);
1126     TwoDPlotFile    = opt2fn_null("-2d", NFILE, fnm);
1127     ThreeDPlotFile  = opt2fn_null("-3d", NFILE, fnm);
1128     FilterFile      = opt2fn_null("-filt", NFILE, fnm);
1129     ExtremeFile     = opt2fn_null("-extr", NFILE, fnm);
1130     OverlapFile     = opt2fn_null("-over", NFILE, fnm);
1131     InpMatFile      = ftp2fn_null(efXPM, NFILE, fnm);
1132
1133     bTop   = fn2bTPX(topfile);
1134     bProj  = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1135         || FilterFile || ExtremeFile;
1136     bFirstLastSet  =
1137         opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1138     bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1139         OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1140     bVec2  = Vec2File || OverlapFile || InpMatFile;
1141     bM     = RmsfFile || bProj;
1142     bTraj  = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1143         || TwoDPlotFile || ThreeDPlotFile;
1144     bIndex = bM || bProj;
1145     bTPS   = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1146         FilterFile  || (bIndex && indexfile);
1147     bCompare = Vec2File || Eig2File;
1148     bPDB3D   = fn2ftp(ThreeDPlotFile) == efPDB;
1149
1150     read_eigenvectors(VecFile, &natoms, &bFit1,
1151                       &xref1, &bDMR1, &xav1, &bDMA1,
1152                       &nvec1, &eignr1, &eigvec1, &eigval1);
1153     neig1 = DIM*natoms;
1154
1155     /* Overwrite eigenvalues from separate files if the user provides them */
1156     if (EigFile != NULL)
1157     {
1158         int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1159         if (neig_tmp != neig1)
1160         {
1161             fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1162         }
1163         neig1 = neig_tmp;
1164         srenew(eigval1, neig1);
1165         for (j = 0; j < neig1; j++)
1166         {
1167             real tmp = eigval1[j];
1168             eigval1[j] = xvgdata[1][j];
1169             if (debug && (eigval1[j] != tmp))
1170             {
1171                 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1172                         j, tmp, eigval1[j]);
1173             }
1174         }
1175         for (j = 0; j < i; j++)
1176         {
1177             sfree(xvgdata[j]);
1178         }
1179         sfree(xvgdata);
1180         fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1181     }
1182
1183     if (bEntropy)
1184     {
1185         if (bDMA1)
1186         {
1187             gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1188         }
1189         calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1190         calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1191     }
1192
1193     if (bVec2)
1194     {
1195         if (!Vec2File)
1196         {
1197             gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1198         }
1199         read_eigenvectors(Vec2File, &neig2, &bFit2,
1200                           &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1201
1202         neig2 = DIM*neig2;
1203         if (neig2 != neig1)
1204         {
1205             gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1206         }
1207     }
1208
1209     if (Eig2File != NULL)
1210     {
1211         neig2 = read_xvg(Eig2File, &xvgdata, &i);
1212         srenew(eigval2, neig2);
1213         for (j = 0; j < neig2; j++)
1214         {
1215             eigval2[j] = xvgdata[1][j];
1216         }
1217         for (j = 0; j < i; j++)
1218         {
1219             sfree(xvgdata[j]);
1220         }
1221         sfree(xvgdata);
1222         fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1223     }
1224
1225
1226     if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1227     {
1228         bM = FALSE;
1229     }
1230     if ((xref1 == NULL) && (bM || bTraj))
1231     {
1232         bTPS = TRUE;
1233     }
1234
1235     xtop  = NULL;
1236     nfit  = 0;
1237     ifit  = NULL;
1238     w_rls = NULL;
1239
1240     if (!bTPS)
1241     {
1242         bTop = FALSE;
1243     }
1244     else
1245     {
1246         bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1247                              title, &top, &ePBC, &xtop, NULL, topbox, bM);
1248         atoms = &top.atoms;
1249         gpbc  = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1250         gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1251         /* Fitting is only required for the projection */
1252         if (bProj && bFit1)
1253         {
1254             if (xref1 == NULL)
1255             {
1256                 printf("\nNote: the structure in %s should be the same\n"
1257                        "      as the one used for the fit in g_covar\n", topfile);
1258             }
1259             printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1260             get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1261
1262             snew(w_rls, atoms->nr);
1263             for (i = 0; (i < nfit); i++)
1264             {
1265                 if (bDMR1)
1266                 {
1267                     w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1268                 }
1269                 else
1270                 {
1271                     w_rls[ifit[i]] = 1.0;
1272                 }
1273             }
1274
1275             snew(xrefp, atoms->nr);
1276             if (xref1 != NULL)
1277             {
1278                 for (i = 0; (i < nfit); i++)
1279                 {
1280                     copy_rvec(xref1[i], xrefp[ifit[i]]);
1281                 }
1282             }
1283             else
1284             {
1285                 /* The top coordinates are the fitting reference */
1286                 for (i = 0; (i < nfit); i++)
1287                 {
1288                     copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1289                 }
1290                 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1291             }
1292         }
1293         gmx_rmpbc_done(gpbc);
1294     }
1295
1296     if (bIndex)
1297     {
1298         printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1299         get_index(atoms, indexfile, 1, &i, &index, &grpname);
1300         if (i != natoms)
1301         {
1302             gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1303         }
1304         printf("\n");
1305     }
1306
1307     snew(sqrtm, natoms);
1308     if (bM && bDMA1)
1309     {
1310         proj_unit = "u\\S1/2\\Nnm";
1311         for (i = 0; (i < natoms); i++)
1312         {
1313             sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1314         }
1315     }
1316     else
1317     {
1318         proj_unit = "nm";
1319         for (i = 0; (i < natoms); i++)
1320         {
1321             sqrtm[i] = 1.0;
1322         }
1323     }
1324
1325     if (bVec2)
1326     {
1327         t       = 0;
1328         totmass = 0;
1329         for (i = 0; (i < natoms); i++)
1330         {
1331             for (d = 0; (d < DIM); d++)
1332             {
1333                 t       += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1334                 totmass += sqr(sqrtm[i]);
1335             }
1336         }
1337         fprintf(stdout, "RMSD (without fit) between the two average structures:"
1338                 " %.3f (nm)\n\n", sqrt(t/totmass));
1339     }
1340
1341     if (last == -1)
1342     {
1343         last = natoms*DIM;
1344     }
1345     if (first > -1)
1346     {
1347         if (bFirstToLast)
1348         {
1349             /* make an index from first to last */
1350             nout = last-first+1;
1351             snew(iout, nout);
1352             for (i = 0; i < nout; i++)
1353             {
1354                 iout[i] = first-1+i;
1355             }
1356         }
1357         else if (ThreeDPlotFile)
1358         {
1359             /* make an index of first+(0,1,2) and last */
1360             nout = bPDB3D ? 4 : 3;
1361             nout = min(last-first+1, nout);
1362             snew(iout, nout);
1363             iout[0] = first-1;
1364             iout[1] = first;
1365             if (nout > 3)
1366             {
1367                 iout[2] = first+1;
1368             }
1369             iout[nout-1] = last-1;
1370         }
1371         else
1372         {
1373             /* make an index of first and last */
1374             nout = 2;
1375             snew(iout, nout);
1376             iout[0] = first-1;
1377             iout[1] = last-1;
1378         }
1379     }
1380     else
1381     {
1382         printf("Select eigenvectors for output, end your selection with 0\n");
1383         nout = -1;
1384         iout = NULL;
1385
1386         do
1387         {
1388             nout++;
1389             srenew(iout, nout+1);
1390             if (1 != scanf("%d", &iout[nout]))
1391             {
1392                 gmx_fatal(FARGS, "Error reading user input");
1393             }
1394             iout[nout]--;
1395         }
1396         while (iout[nout] >= 0);
1397
1398         printf("\n");
1399     }
1400     /* make an index of the eigenvectors which are present */
1401     snew(outvec, nout);
1402     noutvec = 0;
1403     for (i = 0; i < nout; i++)
1404     {
1405         j = 0;
1406         while ((j < nvec1) && (eignr1[j] != iout[i]))
1407         {
1408             j++;
1409         }
1410         if ((j < nvec1) && (eignr1[j] == iout[i]))
1411         {
1412             outvec[noutvec] = j;
1413             noutvec++;
1414         }
1415     }
1416     fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1417     if (noutvec <= 100)
1418     {
1419         fprintf(stderr, ":");
1420         for (j = 0; j < noutvec; j++)
1421         {
1422             fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1423         }
1424     }
1425     fprintf(stderr, "\n");
1426
1427     if (CompFile)
1428     {
1429         components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1430     }
1431
1432     if (RmsfFile)
1433     {
1434         rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1435              neig1, oenv);
1436     }
1437
1438     if (bProj)
1439     {
1440         project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1441                 bTop ? &top : NULL, ePBC, topbox,
1442                 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1443                 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1444                 bFit1, xrefp, nfit, ifit, w_rls,
1445                 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1446                 oenv);
1447     }
1448
1449     if (OverlapFile)
1450     {
1451         overlap(OverlapFile, natoms,
1452                 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1453     }
1454
1455     if (InpMatFile)
1456     {
1457         inprod_matrix(InpMatFile, natoms,
1458                       nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1459                       bFirstLastSet, noutvec, outvec);
1460     }
1461
1462     if (bCompare)
1463     {
1464         compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1465     }
1466
1467
1468     if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1469         !bCompare && !bEntropy)
1470     {
1471         fprintf(stderr, "\nIf you want some output,"
1472                 " set one (or two or ...) of the output file options\n");
1473     }
1474
1475
1476     view_all(oenv, NFILE, fnm);
1477
1478     return 0;
1479 }