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