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