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