src/gromacs/gmxana/gmx_anaeig.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include <math.h>
  40 #include <stdlib.h>
  41 #include <string.h>
  42
  43 #include "gromacs/commandline/pargs.h"
  44 #include "gromacs/fileio/confio.h"
  45 #include "gromacs/fileio/matio.h"
  46 #include "gromacs/fileio/pdbio.h"
  47 #include "gromacs/fileio/tpxio.h"
  48 #include "gromacs/fileio/trxio.h"
  49 #include "gromacs/fileio/xvgr.h"
  50 #include "gromacs/gmxana/eigio.h"
  51 #include "gromacs/gmxana/gmx_ana.h"
  52 #include "gromacs/legacyheaders/macros.h"
  53 #include "gromacs/legacyheaders/txtdump.h"
  54 #include "gromacs/legacyheaders/typedefs.h"
  55 #include "gromacs/legacyheaders/viewit.h"
  56 #include "gromacs/math/do_fit.h"
  57 #include "gromacs/math/units.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/pbcutil/rmpbc.h"
  60 #include "gromacs/topology/index.h"
  61 #include "gromacs/utility/fatalerror.h"
  62 #include "gromacs/utility/futil.h"
  63 #include "gromacs/utility/smalloc.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/gmx_expm1(hwkT) - gmx_log1p(-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 = gmx_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 && fabs(x[i]) < 1e-5)
 259                 {
 260                     fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 261                 }
 262                 fprintf(out, "%10.4f %10.5f\n",
 263                         x[i]*scale_x, y ? y[g][i] : sy[g][s][i]);
 264             }
 265             fprintf(out, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 266         }
 267     }
 268     gmx_ffclose(out);
 269 }
 270
 271 static void
 272 compare(int natoms, int n1, rvec **eigvec1, int n2, rvec **eigvec2,
 273         real *eigval1, int neig1, real *eigval2, int neig2)
 274 {
 275     int    n, nrow;
 276     int    i, j, k;
 277     double sum1, sum2, trace1, trace2, sab, samsb2, tmp, ip;
 278
 279     n = min(n1, n2);
 280
 281     n = min(n, min(neig1, neig2));
 282     fprintf(stdout, "Will compare the covariance matrices using %d dimensions\n", n);
 283
 284     sum1 = 0;
 285     for (i = 0; i < n; i++)
 286     {
 287         if (eigval1[i] < 0)
 288         {
 289             eigval1[i] = 0;
 290         }
 291         sum1      += eigval1[i];
 292         eigval1[i] = sqrt(eigval1[i]);
 293     }
 294     trace1 = sum1;
 295     for (i = n; i < neig1; i++)
 296     {
 297         trace1 += eigval1[i];
 298     }
 299     sum2 = 0;
 300     for (i = 0; i < n; i++)
 301     {
 302         if (eigval2[i] < 0)
 303         {
 304             eigval2[i] = 0;
 305         }
 306         sum2      += eigval2[i];
 307         eigval2[i] = sqrt(eigval2[i]);
 308     }
 309     trace2 = sum2;
 310     for (i = n; i < neig2; i++)
 311     {
 312         trace2 += eigval2[i];
 313     }
 314
 315     fprintf(stdout, "Trace of the two matrices: %g and %g\n", sum1, sum2);
 316     if (neig1 != n || neig2 != n)
 317     {
 318         fprintf(stdout, "this is %d%% and %d%% of the total trace\n",
 319                 (int)(100*sum1/trace1+0.5), (int)(100*sum2/trace2+0.5));
 320     }
 321     fprintf(stdout, "Square root of the traces: %g and %g\n",
 322             sqrt(sum1), sqrt(sum2));
 323
 324     sab = 0;
 325     for (i = 0; i < n; i++)
 326     {
 327         tmp = 0;
 328         for (j = 0; j < n; j++)
 329         {
 330             ip = 0;
 331             for (k = 0; k < natoms; k++)
 332             {
 333                 ip += iprod(eigvec1[i][k], eigvec2[j][k]);
 334             }
 335             tmp += eigval2[j]*ip*ip;
 336         }
 337         sab += eigval1[i]*tmp;
 338     }
 339
 340     samsb2 = sum1+sum2-2*sab;
 341     if (samsb2 < 0)
 342     {
 343         samsb2 = 0;
 344     }
 345
 346     fprintf(stdout, "The overlap of the covariance matrices:\n");
 347     fprintf(stdout, "  normalized:  %.3f\n", 1-sqrt(samsb2/(sum1+sum2)));
 348     tmp = 1-sab/sqrt(sum1*sum2);
 349     if (tmp < 0)
 350     {
 351         tmp = 0;
 352     }
 353     fprintf(stdout, "       shape:  %.3f\n", 1-sqrt(tmp));
 354 }
 355
 356
 357 static void inprod_matrix(const char *matfile, int natoms,
 358                           int nvec1, int *eignr1, rvec **eigvec1,
 359                           int nvec2, int *eignr2, rvec **eigvec2,
 360                           gmx_bool bSelect, int noutvec, int *outvec)
 361 {
 362     FILE   *out;
 363     real  **mat;
 364     int     i, x1, y1, x, y, nlevels;
 365     int     nx, ny;
 366     real    inp, *t_x, *t_y, max;
 367     t_rgb   rlo, rhi;
 368
 369     snew(t_y, nvec2);
 370     if (bSelect)
 371     {
 372         nx = noutvec;
 373         ny = 0;
 374         for (y1 = 0; y1 < nx; y1++)
 375         {
 376             if (outvec[y1] < nvec2)
 377             {
 378                 t_y[ny] = eignr2[outvec[y1]]+1;
 379                 ny++;
 380             }
 381         }
 382     }
 383     else
 384     {
 385         nx = nvec1;
 386         ny = nvec2;
 387         for (y = 0; y < ny; y++)
 388         {
 389             t_y[y] = eignr2[y]+1;
 390         }
 391     }
 392
 393     fprintf(stderr, "Calculating inner-product matrix of %dx%d eigenvectors\n",
 394             nx, nvec2);
 395
 396     snew(mat, nx);
 397     snew(t_x, nx);
 398     max = 0;
 399     for (x1 = 0; x1 < nx; x1++)
 400     {
 401         snew(mat[x1], ny);
 402         if (bSelect)
 403         {
 404             x = outvec[x1];
 405         }
 406         else
 407         {
 408             x = x1;
 409         }
 410         t_x[x1] = eignr1[x]+1;
 411         fprintf(stderr, " %d", eignr1[x]+1);
 412         for (y1 = 0; y1 < ny; y1++)
 413         {
 414             inp = 0;
 415             if (bSelect)
 416             {
 417                 while (outvec[y1] >= nvec2)
 418                 {
 419                     y1++;
 420                 }
 421                 y = outvec[y1];
 422             }
 423             else
 424             {
 425                 y = y1;
 426             }
 427             for (i = 0; i < natoms; i++)
 428             {
 429                 inp += iprod(eigvec1[x][i], eigvec2[y][i]);
 430             }
 431             mat[x1][y1] = fabs(inp);
 432             if (mat[x1][y1] > max)
 433             {
 434                 max = mat[x1][y1];
 435             }
 436         }
 437     }
 438     fprintf(stderr, "\n");
 439     rlo.r   = 1; rlo.g = 1; rlo.b = 1;
 440     rhi.r   = 0; rhi.g = 0; rhi.b = 0;
 441     nlevels = 41;
 442     out     = gmx_ffopen(matfile, "w");
 443     write_xpm(out, 0, "Eigenvector inner-products", "in.prod.", "run 1", "run 2",
 444               nx, ny, t_x, t_y, mat, 0.0, max, rlo, rhi, &nlevels);
 445     gmx_ffclose(out);
 446 }
 447
 448 static void overlap(const char *outfile, int natoms,
 449                     rvec **eigvec1,
 450                     int nvec2, int *eignr2, rvec **eigvec2,
 451                     int noutvec, int *outvec,
 452                     const output_env_t oenv)
 453 {
 454     FILE *out;
 455     int   i, v, vec, x;
 456     real  overlap, inp;
 457
 458     fprintf(stderr, "Calculating overlap between eigenvectors of set 2 with eigenvectors\n");
 459     for (i = 0; i < noutvec; i++)
 460     {
 461         fprintf(stderr, "%d ", outvec[i]+1);
 462     }
 463     fprintf(stderr, "\n");
 464
 465     out = xvgropen(outfile, "Subspace overlap",
 466                    "Eigenvectors of trajectory 2", "Overlap", oenv);
 467     if (output_env_get_print_xvgr_codes(oenv))
 468     {
 469         fprintf(out, "@ subtitle \"using %d eigenvectors of trajectory 1\"\n", noutvec);
 470     }
 471     overlap = 0;
 472     for (x = 0; x < nvec2; x++)
 473     {
 474         for (v = 0; v < noutvec; v++)
 475         {
 476             vec = outvec[v];
 477             inp = 0;
 478             for (i = 0; i < natoms; i++)
 479             {
 480                 inp += iprod(eigvec1[vec][i], eigvec2[x][i]);
 481             }
 482             overlap += sqr(inp);
 483         }
 484         fprintf(out, "%5d  %5.3f\n", eignr2[x]+1, overlap/noutvec);
 485     }
 486
 487     xvgrclose(out);
 488 }
 489
 490 static void project(const char *trajfile, t_topology *top, int ePBC, matrix topbox,
 491                     const char *projfile, const char *twodplotfile,
 492                     const char *threedplotfile, const char *filterfile, int skip,
 493                     const char *extremefile, gmx_bool bExtrAll, real extreme,
 494                     int nextr, t_atoms *atoms, int natoms, atom_id *index,
 495                     gmx_bool bFit, rvec *xref, int nfit, atom_id *ifit, real *w_rls,
 496                     real *sqrtm, rvec *xav,
 497                     int *eignr, rvec **eigvec,
 498                     int noutvec, int *outvec, gmx_bool bSplit,
 499                     const output_env_t oenv)
 500 {
 501     FILE        *xvgrout = NULL;
 502     int          nat, i, j, d, v, vec, nfr, nframes = 0, snew_size, frame;
 503     t_trxstatus *out = NULL;
 504     t_trxstatus *status;
 505     int          noutvec_extr, imin, imax;
 506     real        *pmin, *pmax;
 507     atom_id     *all_at;
 508     matrix       box;
 509     rvec        *xread, *x;
 510     real         t, inp, **inprod = NULL, min = 0, max = 0;
 511     char         str[STRLEN], str2[STRLEN], **ylabel, *c;
 512     real         fact;
 513     gmx_rmpbc_t  gpbc = NULL;
 514
 515     snew(x, natoms);
 516
 517     if (bExtrAll)
 518     {
 519         noutvec_extr = noutvec;
 520     }
 521     else
 522     {
 523         noutvec_extr = 1;
 524     }
 525
 526
 527     if (trajfile)
 528     {
 529         snew(inprod, noutvec+1);
 530
 531         if (filterfile)
 532         {
 533             fprintf(stderr, "Writing a filtered trajectory to %s using eigenvectors\n",
 534                     filterfile);
 535             for (i = 0; i < noutvec; i++)
 536             {
 537                 fprintf(stderr, "%d ", outvec[i]+1);
 538             }
 539             fprintf(stderr, "\n");
 540             out = open_trx(filterfile, "w");
 541         }
 542         snew_size = 0;
 543         nfr       = 0;
 544         nframes   = 0;
 545         nat       = read_first_x(oenv, &status, trajfile, &t, &xread, box);
 546         if (nat > atoms->nr)
 547         {
 548             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);
 549         }
 550         snew(all_at, nat);
 551
 552         if (top)
 553         {
 554             gpbc = gmx_rmpbc_init(&top->idef, ePBC, nat);
 555         }
 556
 557         for (i = 0; i < nat; i++)
 558         {
 559             all_at[i] = i;
 560         }
 561         do
 562         {
 563             if (nfr % skip == 0)
 564             {
 565                 if (top)
 566                 {
 567                     gmx_rmpbc(gpbc, nat, box, xread);
 568                 }
 569                 if (nframes >= snew_size)
 570                 {
 571                     snew_size += 100;
 572                     for (i = 0; i < noutvec+1; i++)
 573                     {
 574                         srenew(inprod[i], snew_size);
 575                     }
 576                 }
 577                 inprod[noutvec][nframes] = t;
 578                 /* calculate x: a fitted struture of the selected atoms */
 579                 if (bFit)
 580                 {
 581                     reset_x(nfit, ifit, nat, NULL, xread, w_rls);
 582                     do_fit(nat, w_rls, xref, xread);
 583                 }
 584                 for (i = 0; i < natoms; i++)
 585                 {
 586                     copy_rvec(xread[index[i]], x[i]);
 587                 }
 588
 589                 for (v = 0; v < noutvec; v++)
 590                 {
 591                     vec = outvec[v];
 592                     /* calculate (mass-weighted) projection */
 593                     inp = 0;
 594                     for (i = 0; i < natoms; i++)
 595                     {
 596                         inp += (eigvec[vec][i][0]*(x[i][0]-xav[i][0])+
 597                                 eigvec[vec][i][1]*(x[i][1]-xav[i][1])+
 598                                 eigvec[vec][i][2]*(x[i][2]-xav[i][2]))*sqrtm[i];
 599                     }
 600                     inprod[v][nframes] = inp;
 601                 }
 602                 if (filterfile)
 603                 {
 604                     for (i = 0; i < natoms; i++)
 605                     {
 606                         for (d = 0; d < DIM; d++)
 607                         {
 608                             /* misuse xread for output */
 609                             xread[index[i]][d] = xav[i][d];
 610                             for (v = 0; v < noutvec; v++)
 611                             {
 612                                 xread[index[i]][d] +=
 613                                     inprod[v][nframes]*eigvec[outvec[v]][i][d]/sqrtm[i];
 614                             }
 615                         }
 616                     }
 617                     write_trx(out, natoms, index, atoms, 0, t, box, xread, NULL, NULL);
 618                 }
 619                 nframes++;
 620             }
 621             nfr++;
 622         }
 623         while (read_next_x(oenv, status, &t, xread, box));
 624         close_trx(status);
 625         sfree(x);
 626         if (filterfile)
 627         {
 628             close_trx(out);
 629         }
 630     }
 631     else
 632     {
 633         snew(xread, atoms->nr);
 634     }
 635
 636     if (top)
 637     {
 638         gmx_rmpbc_done(gpbc);
 639     }
 640
 641
 642     if (projfile)
 643     {
 644         snew(ylabel, noutvec);
 645         for (v = 0; v < noutvec; v++)
 646         {
 647             sprintf(str, "vec %d", eignr[outvec[v]]+1);
 648             ylabel[v] = gmx_strdup(str);
 649         }
 650         sprintf(str, "projection on eigenvectors (%s)", proj_unit);
 651         write_xvgr_graphs(projfile, noutvec, 1, str, NULL, output_env_get_xvgr_tlabel(oenv),
 652                           (const char **)ylabel,
 653                           nframes, inprod[noutvec], inprod, NULL,
 654                           output_env_get_time_factor(oenv), FALSE, bSplit, oenv);
 655     }
 656
 657     if (twodplotfile)
 658     {
 659         sprintf(str, "projection on eigenvector %d (%s)",
 660                 eignr[outvec[0]]+1, proj_unit);
 661         sprintf(str2, "projection on eigenvector %d (%s)",
 662                 eignr[outvec[noutvec-1]]+1, proj_unit);
 663         xvgrout = xvgropen(twodplotfile, "2D projection of trajectory", str, str2,
 664                            oenv);
 665         for (i = 0; i < nframes; i++)
 666         {
 667             if (bSplit && i > 0 && fabs(inprod[noutvec][i]) < 1e-5)
 668             {
 669                 fprintf(xvgrout, "%s\n", output_env_get_print_xvgr_codes(oenv) ? "&" : "");
 670             }
 671             fprintf(xvgrout, "%10.5f %10.5f\n", inprod[0][i], inprod[noutvec-1][i]);
 672         }
 673         xvgrclose(xvgrout);
 674     }
 675
 676     if (threedplotfile)
 677     {
 678         t_atoms     atoms;
 679         rvec       *x;
 680         real       *b = NULL;
 681         matrix      box;
 682         char       *resnm, *atnm;
 683         gmx_bool    bPDB, b4D;
 684         FILE       *out;
 685
 686         if (noutvec < 3)
 687         {
 688             gmx_fatal(FARGS, "You have selected less than 3 eigenvectors");
 689         }
 690
 691         /* initialize */
 692         bPDB = fn2ftp(threedplotfile) == efPDB;
 693         clear_mat(box);
 694         box[XX][XX] = box[YY][YY] = box[ZZ][ZZ] = 1;
 695
 696         b4D = bPDB && (noutvec >= 4);
 697         if (b4D)
 698         {
 699             fprintf(stderr, "You have selected four or more eigenvectors:\n"
 700                     "fourth eigenvector will be plotted "
 701                     "in bfactor field of pdb file\n");
 702             sprintf(str, "4D proj. of traj. on eigenv. %d, %d, %d and %d",
 703                     eignr[outvec[0]]+1, eignr[outvec[1]]+1,
 704                     eignr[outvec[2]]+1, eignr[outvec[3]]+1);
 705         }
 706         else
 707         {
 708             sprintf(str, "3D proj. of traj. on eigenv. %d, %d and %d",
 709                     eignr[outvec[0]]+1, eignr[outvec[1]]+1, eignr[outvec[2]]+1);
 710         }
 711         init_t_atoms(&atoms, nframes, FALSE);
 712         snew(x, nframes);
 713         snew(b, nframes);
 714         atnm  = gmx_strdup("C");
 715         resnm = gmx_strdup("PRJ");
 716
 717         if (nframes > 10000)
 718         {
 719             fact = 10000.0/nframes;
 720         }
 721         else
 722         {
 723             fact = 1.0;
 724         }
 725
 726         for (i = 0; i < nframes; i++)
 727         {
 728             atoms.atomname[i]     = &atnm;
 729             atoms.atom[i].resind  = i;
 730             atoms.resinfo[i].name = &resnm;
 731             atoms.resinfo[i].nr   = ceil(i*fact);
 732             atoms.resinfo[i].ic   = ' ';
 733             x[i][XX]              = inprod[0][i];
 734             x[i][YY]              = inprod[1][i];
 735             x[i][ZZ]              = inprod[2][i];
 736             if (b4D)
 737             {
 738                 b[i]  = inprod[3][i];
 739             }
 740         }
 741         if ( ( b4D || bSplit ) && bPDB)
 742         {
 743             out = gmx_ffopen(threedplotfile, "w");
 744             fprintf(out, "HEADER    %s\n", str);
 745             if (b4D)
 746             {
 747                 fprintf(out, "REMARK    %s\n", "fourth dimension plotted as B-factor");
 748             }
 749             j = 0;
 750             for (i = 0; i < atoms.nr; i++)
 751             {
 752                 if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
 753                 {
 754                     fprintf(out, "TER\n");
 755                     j = 0;
 756                 }
 757                 gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
 758                                          10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
 759                 if (j > 0)
 760                 {
 761                     fprintf(out, "CONECT%5d%5d\n", i, i+1);
 762                 }
 763                 j++;
 764             }
 765             fprintf(out, "TER\n");
 766             gmx_ffclose(out);
 767         }
 768         else
 769         {
 770             write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
 771         }
 772         free_t_atoms(&atoms, FALSE);
 773     }
 774
 775     if (extremefile)
 776     {
 777         snew(pmin, noutvec_extr);
 778         snew(pmax, noutvec_extr);
 779         if (extreme == 0)
 780         {
 781             fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
 782             fprintf(stderr,
 783                     "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
 784             imin = 0;
 785             imax = 0;
 786             for (v = 0; v < noutvec_extr; v++)
 787             {
 788                 for (i = 0; i < nframes; i++)
 789                 {
 790                     if (inprod[v][i] < inprod[v][imin])
 791                     {
 792                         imin = i;
 793                     }
 794                     if (inprod[v][i] > inprod[v][imax])
 795                     {
 796                         imax = i;
 797                     }
 798                 }
 799                 pmin[v] = inprod[v][imin];
 800                 pmax[v] = inprod[v][imax];
 801                 fprintf(stderr, "%7d     %10.6f %10d %10.6f %10d\n",
 802                         eignr[outvec[v]]+1,
 803                         pmin[v], imin, pmax[v], imax);
 804             }
 805         }
 806         else
 807         {
 808             pmin[0] = -extreme;
 809             pmax[0] =  extreme;
 810         }
 811         /* build format string for filename: */
 812         strcpy(str, extremefile);  /* copy filename */
 813         c = strrchr(str, '.');     /* find where extention begins */
 814         strcpy(str2, c);           /* get extention */
 815         sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
 816         for (v = 0; v < noutvec_extr; v++)
 817         {
 818             /* make filename using format string */
 819             if (noutvec_extr == 1)
 820             {
 821                 strcpy(str2, extremefile);
 822             }
 823             else
 824             {
 825                 sprintf(str2, str, eignr[outvec[v]]+1);
 826             }
 827             fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
 828                     nextr, outvec[v]+1, str2);
 829             out = open_trx(str2, "w");
 830             for (frame = 0; frame < nextr; frame++)
 831             {
 832                 if ((extreme == 0) && (nextr <= 3))
 833                 {
 834                     for (i = 0; i < natoms; i++)
 835                     {
 836                         atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
 837                     }
 838                 }
 839                 for (i = 0; i < natoms; i++)
 840                 {
 841                     for (d = 0; d < DIM; d++)
 842                     {
 843                         xread[index[i]][d] =
 844                             (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
 845                              *eigvec[outvec[v]][i][d]/sqrtm[i]);
 846                     }
 847                 }
 848                 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
 849             }
 850             close_trx(out);
 851         }
 852         sfree(pmin);
 853         sfree(pmax);
 854     }
 855     fprintf(stderr, "\n");
 856 }
 857
 858 static void components(const char *outfile, int natoms,
 859                        int *eignr, rvec **eigvec,
 860                        int noutvec, int *outvec,
 861                        const output_env_t oenv)
 862 {
 863     int   g, s, v, i;
 864     real *x, ***y;
 865     char  str[STRLEN], **ylabel;
 866
 867     fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
 868
 869     snew(ylabel, noutvec);
 870     snew(y, noutvec);
 871     snew(x, natoms);
 872     for (i = 0; i < natoms; i++)
 873     {
 874         x[i] = i+1;
 875     }
 876     for (g = 0; g < noutvec; g++)
 877     {
 878         v = outvec[g];
 879         sprintf(str, "vec %d", eignr[v]+1);
 880         ylabel[g] = gmx_strdup(str);
 881         snew(y[g], 4);
 882         for (s = 0; s < 4; s++)
 883         {
 884             snew(y[g][s], natoms);
 885         }
 886         for (i = 0; i < natoms; i++)
 887         {
 888             y[g][0][i] = norm(eigvec[v][i]);
 889             for (s = 0; s < 3; s++)
 890             {
 891                 y[g][s+1][i] = eigvec[v][i][s];
 892             }
 893         }
 894     }
 895     write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
 896                       "black: total, red: x, green: y, blue: z",
 897                       "Atom number", (const char **)ylabel,
 898                       natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
 899     fprintf(stderr, "\n");
 900 }
 901
 902 static void rmsf(const char *outfile, int natoms, real *sqrtm,
 903                  int *eignr, rvec **eigvec,
 904                  int noutvec, int *outvec,
 905                  real *eigval, int neig,
 906                  const output_env_t oenv)
 907 {
 908     int    nrow, g, v, i;
 909     real  *x, **y;
 910     char   str[STRLEN], **ylabel;
 911
 912     for (i = 0; i < neig; i++)
 913     {
 914         if (eigval[i] < 0)
 915         {
 916             eigval[i] = 0;
 917         }
 918     }
 919
 920     fprintf(stderr, "Writing rmsf to %s\n", outfile);
 921
 922     snew(ylabel, noutvec);
 923     snew(y, noutvec);
 924     snew(x, natoms);
 925     for (i = 0; i < natoms; i++)
 926     {
 927         x[i] = i+1;
 928     }
 929     for (g = 0; g < noutvec; g++)
 930     {
 931         v = outvec[g];
 932         if (eignr[v] >= neig)
 933         {
 934             gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
 935         }
 936         sprintf(str, "vec %d", eignr[v]+1);
 937         ylabel[g] = gmx_strdup(str);
 938         snew(y[g], natoms);
 939         for (i = 0; i < natoms; i++)
 940         {
 941             y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
 942         }
 943     }
 944     write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
 945                       "Atom number", (const char **)ylabel,
 946                       natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
 947     fprintf(stderr, "\n");
 948 }
 949
 950 int gmx_anaeig(int argc, char *argv[])
 951 {
 952     static const char *desc[] = {
 953         "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
 954         "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
 955         "([gmx-nmeig]).[PAR]",
 956
 957         "When a trajectory is projected on eigenvectors, all structures are",
 958         "fitted to the structure in the eigenvector file, if present, otherwise",
 959         "to the structure in the structure file. When no run input file is",
 960         "supplied, periodicity will not be taken into account. Most analyses",
 961         "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
 962         "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
 963
 964         "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
 965         "[TT]-first[tt] to [TT]-last[tt].[PAR]",
 966
 967         "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
 968         "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
 969
 970         "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
 971         "[TT]-first[tt] to [TT]-last[tt].",
 972         "The projections of a trajectory on the eigenvectors of its",
 973         "covariance matrix are called principal components (pc's).",
 974         "It is often useful to check the cosine content of the pc's,",
 975         "since the pc's of random diffusion are cosines with the number",
 976         "of periods equal to half the pc index.",
 977         "The cosine content of the pc's can be calculated with the program",
 978         "[gmx-analyze].[PAR]",
 979
 980         "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
 981         "[TT]-first[tt] and [TT]-last[tt].[PAR]",
 982
 983         "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
 984         "three selected eigenvectors.[PAR]",
 985
 986         "[TT]-filt[tt]: filter the trajectory to show only the motion along",
 987         "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
 988
 989         "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
 990         "on the average structure and interpolate [TT]-nframes[tt] frames",
 991         "between them, or set your own extremes with [TT]-max[tt]. The",
 992         "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
 993         "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
 994         "will be written to separate files. Chain identifiers will be added",
 995         "when writing a [REF].pdb[ref] file with two or three structures (you",
 996         "can use [TT]rasmol -nmrpdb[tt] to view such a [REF].pdb[ref] file).[PAR]",
 997
 998         "  Overlap calculations between covariance analysis:[BR]",
 999         "  [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1000
1001         "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1002         "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1003         "in file [TT]-v[tt].[PAR]",
1004
1005         "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1006         "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1007         "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1008         "have been set explicitly.[PAR]",
1009
1010         "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1011         "a single number for the overlap between the covariance matrices is",
1012         "generated. The formulas are::",
1013         "",
1014         "         difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))",
1015         " normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))",
1016         "      shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))",
1017         "",
1018         "where M1 and M2 are the two covariance matrices and tr is the trace",
1019         "of a matrix. The numbers are proportional to the overlap of the square",
1020         "root of the fluctuations. The normalized overlap is the most useful",
1021         "number, it is 1 for identical matrices and 0 when the sampled",
1022         "subspaces are orthogonal.[PAR]",
1023         "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1024         "computed based on the Quasiharmonic approach and based on",
1025         "Schlitter's formula."
1026     };
1027     static int         first  = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1028     static real        max    = 0.0, temp = 298.15;
1029     static gmx_bool    bSplit = FALSE, bEntropy = FALSE;
1030     t_pargs            pa[]   = {
1031         { "-first", FALSE, etINT, {&first},
1032           "First eigenvector for analysis (-1 is select)" },
1033         { "-last",  FALSE, etINT, {&last},
1034           "Last eigenvector for analysis (-1 is till the last)" },
1035         { "-skip",  FALSE, etINT, {&skip},
1036           "Only analyse every nr-th frame" },
1037         { "-max",  FALSE, etREAL, {&max},
1038           "Maximum for projection of the eigenvector on the average structure, "
1039           "max=0 gives the extremes" },
1040         { "-nframes",  FALSE, etINT, {&nextr},
1041           "Number of frames for the extremes output" },
1042         { "-split",   FALSE, etBOOL, {&bSplit},
1043           "Split eigenvector projections where time is zero" },
1044         { "-entropy", FALSE, etBOOL, {&bEntropy},
1045           "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1046         { "-temp",    FALSE, etREAL, {&temp},
1047           "Temperature for entropy calculations" },
1048         { "-nevskip", FALSE, etINT, {&nskip},
1049           "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." }
1050     };
1051 #define NPA asize(pa)
1052
1053     FILE          *out;
1054     int            status, trjout;
1055     t_topology     top;
1056     int            ePBC  = -1;
1057     t_atoms       *atoms = NULL;
1058     rvec          *xtop, *xref1, *xref2, *xrefp = NULL;
1059     gmx_bool       bDMR1, bDMA1, bDMR2, bDMA2;
1060     int            nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1061     rvec          *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1062     matrix         topbox;
1063     real           xid, totmass, *sqrtm, *w_rls, t, lambda;
1064     int            natoms, step;
1065     char          *grpname;
1066     const char    *indexfile;
1067     char           title[STRLEN];
1068     int            i, j, d;
1069     int            nout, *iout, noutvec, *outvec, nfit;
1070     atom_id       *index, *ifit;
1071     const char    *VecFile, *Vec2File, *topfile;
1072     const char    *EigFile, *Eig2File;
1073     const char    *CompFile, *RmsfFile, *ProjOnVecFile;
1074     const char    *TwoDPlotFile, *ThreeDPlotFile;
1075     const char    *FilterFile, *ExtremeFile;
1076     const char    *OverlapFile, *InpMatFile;
1077     gmx_bool       bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1078     gmx_bool       bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1079     real          *eigval1 = NULL, *eigval2 = NULL;
1080     int            neig1, neig2;
1081     double       **xvgdata;
1082     output_env_t   oenv;
1083     gmx_rmpbc_t    gpbc;
1084
1085     t_filenm       fnm[] = {
1086         { efTRN, "-v",    "eigenvec",    ffREAD  },
1087         { efTRN, "-v2",   "eigenvec2",   ffOPTRD },
1088         { efTRX, "-f",    NULL,          ffOPTRD },
1089         { efTPS, NULL,    NULL,          ffOPTRD },
1090         { efNDX, NULL,    NULL,          ffOPTRD },
1091         { efXVG, "-eig", "eigenval",     ffOPTRD },
1092         { efXVG, "-eig2", "eigenval2",   ffOPTRD },
1093         { efXVG, "-comp", "eigcomp",     ffOPTWR },
1094         { efXVG, "-rmsf", "eigrmsf",     ffOPTWR },
1095         { efXVG, "-proj", "proj",        ffOPTWR },
1096         { efXVG, "-2d",   "2dproj",      ffOPTWR },
1097         { efSTO, "-3d",   "3dproj.pdb",  ffOPTWR },
1098         { efTRX, "-filt", "filtered",    ffOPTWR },
1099         { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1100         { efXVG, "-over", "overlap",     ffOPTWR },
1101         { efXPM, "-inpr", "inprod",      ffOPTWR }
1102     };
1103 #define NFILE asize(fnm)
1104
1105     if (!parse_common_args(&argc, argv,
1106                            PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW,
1107                            NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1108     {
1109         return 0;
1110     }
1111
1112     indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1113
1114     VecFile         = opt2fn("-v", NFILE, fnm);
1115     Vec2File        = opt2fn_null("-v2", NFILE, fnm);
1116     topfile         = ftp2fn(efTPS, NFILE, fnm);
1117     EigFile         = opt2fn_null("-eig", NFILE, fnm);
1118     Eig2File        = opt2fn_null("-eig2", NFILE, fnm);
1119     CompFile        = opt2fn_null("-comp", NFILE, fnm);
1120     RmsfFile        = opt2fn_null("-rmsf", NFILE, fnm);
1121     ProjOnVecFile   = opt2fn_null("-proj", NFILE, fnm);
1122     TwoDPlotFile    = opt2fn_null("-2d", NFILE, fnm);
1123     ThreeDPlotFile  = opt2fn_null("-3d", NFILE, fnm);
1124     FilterFile      = opt2fn_null("-filt", NFILE, fnm);
1125     ExtremeFile     = opt2fn_null("-extr", NFILE, fnm);
1126     OverlapFile     = opt2fn_null("-over", NFILE, fnm);
1127     InpMatFile      = ftp2fn_null(efXPM, NFILE, fnm);
1128
1129     bTop   = fn2bTPX(topfile);
1130     bProj  = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1131         || FilterFile || ExtremeFile;
1132     bFirstLastSet  =
1133         opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1134     bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1135         OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1136     bVec2  = Vec2File || OverlapFile || InpMatFile;
1137     bM     = RmsfFile || bProj;
1138     bTraj  = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1139         || TwoDPlotFile || ThreeDPlotFile;
1140     bIndex = bM || bProj;
1141     bTPS   = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1142         FilterFile  || (bIndex && indexfile);
1143     bCompare = Vec2File || Eig2File;
1144     bPDB3D   = fn2ftp(ThreeDPlotFile) == efPDB;
1145
1146     read_eigenvectors(VecFile, &natoms, &bFit1,
1147                       &xref1, &bDMR1, &xav1, &bDMA1,
1148                       &nvec1, &eignr1, &eigvec1, &eigval1);
1149     neig1 = DIM*natoms;
1150
1151     /* Overwrite eigenvalues from separate files if the user provides them */
1152     if (EigFile != NULL)
1153     {
1154         int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1155         if (neig_tmp != neig1)
1156         {
1157             fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1158         }
1159         neig1 = neig_tmp;
1160         srenew(eigval1, neig1);
1161         for (j = 0; j < neig1; j++)
1162         {
1163             real tmp = eigval1[j];
1164             eigval1[j] = xvgdata[1][j];
1165             if (debug && (eigval1[j] != tmp))
1166             {
1167                 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1168                         j, tmp, eigval1[j]);
1169             }
1170         }
1171         for (j = 0; j < i; j++)
1172         {
1173             sfree(xvgdata[j]);
1174         }
1175         sfree(xvgdata);
1176         fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1177     }
1178
1179     if (bEntropy)
1180     {
1181         if (bDMA1)
1182         {
1183             gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1184         }
1185         calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1186         calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1187     }
1188
1189     if (bVec2)
1190     {
1191         if (!Vec2File)
1192         {
1193             gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1194         }
1195         read_eigenvectors(Vec2File, &neig2, &bFit2,
1196                           &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1197
1198         neig2 = DIM*neig2;
1199         if (neig2 != neig1)
1200         {
1201             gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1202         }
1203     }
1204
1205     if (Eig2File != NULL)
1206     {
1207         neig2 = read_xvg(Eig2File, &xvgdata, &i);
1208         srenew(eigval2, neig2);
1209         for (j = 0; j < neig2; j++)
1210         {
1211             eigval2[j] = xvgdata[1][j];
1212         }
1213         for (j = 0; j < i; j++)
1214         {
1215             sfree(xvgdata[j]);
1216         }
1217         sfree(xvgdata);
1218         fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1219     }
1220
1221
1222     if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1223     {
1224         bM = FALSE;
1225     }
1226     if ((xref1 == NULL) && (bM || bTraj))
1227     {
1228         bTPS = TRUE;
1229     }
1230
1231     xtop  = NULL;
1232     nfit  = 0;
1233     ifit  = NULL;
1234     w_rls = NULL;
1235
1236     if (!bTPS)
1237     {
1238         bTop = FALSE;
1239     }
1240     else
1241     {
1242         bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1243                              title, &top, &ePBC, &xtop, NULL, topbox, bM);
1244         atoms = &top.atoms;
1245         gpbc  = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1246         gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1247         /* Fitting is only required for the projection */
1248         if (bProj && bFit1)
1249         {
1250             if (xref1 == NULL)
1251             {
1252                 printf("\nNote: the structure in %s should be the same\n"
1253                        "      as the one used for the fit in g_covar\n", topfile);
1254             }
1255             printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1256             get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1257
1258             snew(w_rls, atoms->nr);
1259             for (i = 0; (i < nfit); i++)
1260             {
1261                 if (bDMR1)
1262                 {
1263                     w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1264                 }
1265                 else
1266                 {
1267                     w_rls[ifit[i]] = 1.0;
1268                 }
1269             }
1270
1271             snew(xrefp, atoms->nr);
1272             if (xref1 != NULL)
1273             {
1274                 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1275                 if (natoms != nfit)
1276                 {
1277                     gmx_fatal(FARGS, "you selected a group with %d elements instead of %d, your selection does not fit the reference structure in the eigenvector file.", nfit, natoms);
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 }