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 && fabs(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 && fabs(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;
 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             out = gmx_ffopen(threedplotfile, "w");
 750             fprintf(out, "HEADER    %s\n", str);
 751             if (b4D)
 752             {
 753                 fprintf(out, "REMARK    %s\n", "fourth dimension plotted as B-factor");
 754             }
 755             j = 0;
 756             for (i = 0; i < atoms.nr; i++)
 757             {
 758                 if (j > 0 && bSplit && fabs(inprod[noutvec][i]) < 1e-5)
 759                 {
 760                     fprintf(out, "TER\n");
 761                     j = 0;
 762                 }
 763                 gmx_fprintf_pdb_atomline(out, epdbATOM, i+1, "C", ' ', "PRJ", ' ', j+1, ' ',
 764                                          10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 10*b[i], "");
 765                 if (j > 0)
 766                 {
 767                     fprintf(out, "CONECT%5d%5d\n", i, i+1);
 768                 }
 769                 j++;
 770             }
 771             fprintf(out, "TER\n");
 772             gmx_ffclose(out);
 773         }
 774         else
 775         {
 776             write_sto_conf(threedplotfile, str, &atoms, x, NULL, ePBC, box);
 777         }
 778         free_t_atoms(&atoms, FALSE);
 779     }
 780
 781     if (extremefile)
 782     {
 783         snew(pmin, noutvec_extr);
 784         snew(pmax, noutvec_extr);
 785         if (extreme == 0)
 786         {
 787             fprintf(stderr, "%11s %17s %17s\n", "eigenvector", "Minimum", "Maximum");
 788             fprintf(stderr,
 789                     "%11s %10s %10s %10s %10s\n", "", "value", "frame", "value", "frame");
 790             imin = 0;
 791             imax = 0;
 792             for (v = 0; v < noutvec_extr; v++)
 793             {
 794                 for (i = 0; i < nframes; i++)
 795                 {
 796                     if (inprod[v][i] < inprod[v][imin])
 797                     {
 798                         imin = i;
 799                     }
 800                     if (inprod[v][i] > inprod[v][imax])
 801                     {
 802                         imax = i;
 803                     }
 804                 }
 805                 pmin[v] = inprod[v][imin];
 806                 pmax[v] = inprod[v][imax];
 807                 fprintf(stderr, "%7d     %10.6f %10d %10.6f %10d\n",
 808                         eignr[outvec[v]]+1,
 809                         pmin[v], imin, pmax[v], imax);
 810             }
 811         }
 812         else
 813         {
 814             pmin[0] = -extreme;
 815             pmax[0] =  extreme;
 816         }
 817         /* build format string for filename: */
 818         strcpy(str, extremefile);  /* copy filename */
 819         c = strrchr(str, '.');     /* find where extention begins */
 820         strcpy(str2, c);           /* get extention */
 821         sprintf(c, "%%d%s", str2); /* append '%s' and extention to filename */
 822         for (v = 0; v < noutvec_extr; v++)
 823         {
 824             /* make filename using format string */
 825             if (noutvec_extr == 1)
 826             {
 827                 strcpy(str2, extremefile);
 828             }
 829             else
 830             {
 831                 sprintf(str2, str, eignr[outvec[v]]+1);
 832             }
 833             fprintf(stderr, "Writing %d frames along eigenvector %d to %s\n",
 834                     nextr, outvec[v]+1, str2);
 835             out = open_trx(str2, "w");
 836             for (frame = 0; frame < nextr; frame++)
 837             {
 838                 if ((extreme == 0) && (nextr <= 3))
 839                 {
 840                     for (i = 0; i < natoms; i++)
 841                     {
 842                         atoms->resinfo[atoms->atom[index[i]].resind].chainid = 'A' + frame;
 843                     }
 844                 }
 845                 for (i = 0; i < natoms; i++)
 846                 {
 847                     for (d = 0; d < DIM; d++)
 848                     {
 849                         xread[index[i]][d] =
 850                             (xav[i][d] + (pmin[v]*(nextr-frame-1)+pmax[v]*frame)/(nextr-1)
 851                              *eigvec[outvec[v]][i][d]/sqrtm[i]);
 852                     }
 853                 }
 854                 write_trx(out, natoms, index, atoms, 0, frame, topbox, xread, NULL, NULL);
 855             }
 856             close_trx(out);
 857         }
 858         sfree(pmin);
 859         sfree(pmax);
 860     }
 861     fprintf(stderr, "\n");
 862 }
 863
 864 static void components(const char *outfile, int natoms,
 865                        int *eignr, rvec **eigvec,
 866                        int noutvec, int *outvec,
 867                        const output_env_t oenv)
 868 {
 869     int   g, s, v, i;
 870     real *x, ***y;
 871     char  str[STRLEN], **ylabel;
 872
 873     fprintf(stderr, "Writing eigenvector components to %s\n", outfile);
 874
 875     snew(ylabel, noutvec);
 876     snew(y, noutvec);
 877     snew(x, natoms);
 878     for (i = 0; i < natoms; i++)
 879     {
 880         x[i] = i+1;
 881     }
 882     for (g = 0; g < noutvec; g++)
 883     {
 884         v = outvec[g];
 885         sprintf(str, "vec %d", eignr[v]+1);
 886         ylabel[g] = strdup(str);
 887         snew(y[g], 4);
 888         for (s = 0; s < 4; s++)
 889         {
 890             snew(y[g][s], natoms);
 891         }
 892         for (i = 0; i < natoms; i++)
 893         {
 894             y[g][0][i] = norm(eigvec[v][i]);
 895             for (s = 0; s < 3; s++)
 896             {
 897                 y[g][s+1][i] = eigvec[v][i][s];
 898             }
 899         }
 900     }
 901     write_xvgr_graphs(outfile, noutvec, 4, "Eigenvector components",
 902                       "black: total, red: x, green: y, blue: z",
 903                       "Atom number", (const char **)ylabel,
 904                       natoms, x, NULL, y, 1, FALSE, FALSE, oenv);
 905     fprintf(stderr, "\n");
 906 }
 907
 908 static void rmsf(const char *outfile, int natoms, real *sqrtm,
 909                  int *eignr, rvec **eigvec,
 910                  int noutvec, int *outvec,
 911                  real *eigval, int neig,
 912                  const output_env_t oenv)
 913 {
 914     int    nrow, g, v, i;
 915     real  *x, **y;
 916     char   str[STRLEN], **ylabel;
 917
 918     for (i = 0; i < neig; i++)
 919     {
 920         if (eigval[i] < 0)
 921         {
 922             eigval[i] = 0;
 923         }
 924     }
 925
 926     fprintf(stderr, "Writing rmsf to %s\n", outfile);
 927
 928     snew(ylabel, noutvec);
 929     snew(y, noutvec);
 930     snew(x, natoms);
 931     for (i = 0; i < natoms; i++)
 932     {
 933         x[i] = i+1;
 934     }
 935     for (g = 0; g < noutvec; g++)
 936     {
 937         v = outvec[g];
 938         if (eignr[v] >= neig)
 939         {
 940             gmx_fatal(FARGS, "Selected vector %d is larger than the number of eigenvalues (%d)", eignr[v]+1, neig);
 941         }
 942         sprintf(str, "vec %d", eignr[v]+1);
 943         ylabel[g] = strdup(str);
 944         snew(y[g], natoms);
 945         for (i = 0; i < natoms; i++)
 946         {
 947             y[g][i] = sqrt(eigval[eignr[v]]*norm2(eigvec[v][i]))/sqrtm[i];
 948         }
 949     }
 950     write_xvgr_graphs(outfile, noutvec, 1, "RMS fluctuation (nm) ", NULL,
 951                       "Atom number", (const char **)ylabel,
 952                       natoms, x, y, NULL, 1, TRUE, FALSE, oenv);
 953     fprintf(stderr, "\n");
 954 }
 955
 956 int gmx_anaeig(int argc, char *argv[])
 957 {
 958     static const char *desc[] = {
 959         "[THISMODULE] analyzes eigenvectors. The eigenvectors can be of a",
 960         "covariance matrix ([gmx-covar]) or of a Normal Modes analysis",
 961         "([gmx-nmeig]).[PAR]",
 962
 963         "When a trajectory is projected on eigenvectors, all structures are",
 964         "fitted to the structure in the eigenvector file, if present, otherwise",
 965         "to the structure in the structure file. When no run input file is",
 966         "supplied, periodicity will not be taken into account. Most analyses",
 967         "are performed on eigenvectors [TT]-first[tt] to [TT]-last[tt], but when",
 968         "[TT]-first[tt] is set to -1 you will be prompted for a selection.[PAR]",
 969
 970         "[TT]-comp[tt]: plot the vector components per atom of eigenvectors",
 971         "[TT]-first[tt] to [TT]-last[tt].[PAR]",
 972
 973         "[TT]-rmsf[tt]: plot the RMS fluctuation per atom of eigenvectors",
 974         "[TT]-first[tt] to [TT]-last[tt] (requires [TT]-eig[tt]).[PAR]",
 975
 976         "[TT]-proj[tt]: calculate projections of a trajectory on eigenvectors",
 977         "[TT]-first[tt] to [TT]-last[tt].",
 978         "The projections of a trajectory on the eigenvectors of its",
 979         "covariance matrix are called principal components (pc's).",
 980         "It is often useful to check the cosine content of the pc's,",
 981         "since the pc's of random diffusion are cosines with the number",
 982         "of periods equal to half the pc index.",
 983         "The cosine content of the pc's can be calculated with the program",
 984         "[gmx-analyze].[PAR]",
 985
 986         "[TT]-2d[tt]: calculate a 2d projection of a trajectory on eigenvectors",
 987         "[TT]-first[tt] and [TT]-last[tt].[PAR]",
 988
 989         "[TT]-3d[tt]: calculate a 3d projection of a trajectory on the first",
 990         "three selected eigenvectors.[PAR]",
 991
 992         "[TT]-filt[tt]: filter the trajectory to show only the motion along",
 993         "eigenvectors [TT]-first[tt] to [TT]-last[tt].[PAR]",
 994
 995         "[TT]-extr[tt]: calculate the two extreme projections along a trajectory",
 996         "on the average structure and interpolate [TT]-nframes[tt] frames",
 997         "between them, or set your own extremes with [TT]-max[tt]. The",
 998         "eigenvector [TT]-first[tt] will be written unless [TT]-first[tt] and",
 999         "[TT]-last[tt] have been set explicitly, in which case all eigenvectors",
1000         "will be written to separate files. Chain identifiers will be added",
1001         "when writing a [TT].pdb[tt] file with two or three structures (you",
1002         "can use [TT]rasmol -nmrpdb[tt] to view such a [TT].pdb[tt] file).[PAR]",
1003
1004         "  Overlap calculations between covariance analysis:[BR]",
1005         "  [BB]Note:[bb] the analysis should use the same fitting structure[PAR]",
1006
1007         "[TT]-over[tt]: calculate the subspace overlap of the eigenvectors in",
1008         "file [TT]-v2[tt] with eigenvectors [TT]-first[tt] to [TT]-last[tt]",
1009         "in file [TT]-v[tt].[PAR]",
1010
1011         "[TT]-inpr[tt]: calculate a matrix of inner-products between",
1012         "eigenvectors in files [TT]-v[tt] and [TT]-v2[tt]. All eigenvectors",
1013         "of both files will be used unless [TT]-first[tt] and [TT]-last[tt]",
1014         "have been set explicitly.[PAR]",
1015
1016         "When [TT]-v[tt], [TT]-eig[tt], [TT]-v2[tt] and [TT]-eig2[tt] are given,",
1017         "a single number for the overlap between the covariance matrices is",
1018         "generated. The formulas are:[BR]",
1019         "        difference = sqrt(tr((sqrt(M1) - sqrt(M2))^2))[BR]",
1020         "normalized overlap = 1 - difference/sqrt(tr(M1) + tr(M2))[BR]",
1021         "     shape overlap = 1 - sqrt(tr((sqrt(M1/tr(M1)) - sqrt(M2/tr(M2)))^2))[BR]",
1022         "where M1 and M2 are the two covariance matrices and tr is the trace",
1023         "of a matrix. The numbers are proportional to the overlap of the square",
1024         "root of the fluctuations. The normalized overlap is the most useful",
1025         "number, it is 1 for identical matrices and 0 when the sampled",
1026         "subspaces are orthogonal.[PAR]",
1027         "When the [TT]-entropy[tt] flag is given an entropy estimate will be",
1028         "computed based on the Quasiharmonic approach and based on",
1029         "Schlitter's formula."
1030     };
1031     static int         first  = 1, last = -1, skip = 1, nextr = 2, nskip = 6;
1032     static real        max    = 0.0, temp = 298.15;
1033     static gmx_bool    bSplit = FALSE, bEntropy = FALSE;
1034     t_pargs            pa[]   = {
1035         { "-first", FALSE, etINT, {&first},
1036           "First eigenvector for analysis (-1 is select)" },
1037         { "-last",  FALSE, etINT, {&last},
1038           "Last eigenvector for analysis (-1 is till the last)" },
1039         { "-skip",  FALSE, etINT, {&skip},
1040           "Only analyse every nr-th frame" },
1041         { "-max",  FALSE, etREAL, {&max},
1042           "Maximum for projection of the eigenvector on the average structure, "
1043           "max=0 gives the extremes" },
1044         { "-nframes",  FALSE, etINT, {&nextr},
1045           "Number of frames for the extremes output" },
1046         { "-split",   FALSE, etBOOL, {&bSplit},
1047           "Split eigenvector projections where time is zero" },
1048         { "-entropy", FALSE, etBOOL, {&bEntropy},
1049           "Compute entropy according to the Quasiharmonic formula or Schlitter's method." },
1050         { "-temp",    FALSE, etREAL, {&temp},
1051           "Temperature for entropy calculations" },
1052         { "-nevskip", FALSE, etINT, {&nskip},
1053           "Number of eigenvalues to skip when computing the entropy due to the quasi harmonic approximation. When you do a rotational and/or translational fit prior to the covariance analysis, you get 3 or 6 eigenvalues that are very close to zero, and which should not be taken into account when computing the entropy." }
1054     };
1055 #define NPA asize(pa)
1056
1057     FILE          *out;
1058     int            status, trjout;
1059     t_topology     top;
1060     int            ePBC  = -1;
1061     t_atoms       *atoms = NULL;
1062     rvec          *xtop, *xref1, *xref2, *xrefp = NULL;
1063     gmx_bool       bDMR1, bDMA1, bDMR2, bDMA2;
1064     int            nvec1, nvec2, *eignr1 = NULL, *eignr2 = NULL;
1065     rvec          *x, *xread, *xav1, *xav2, **eigvec1 = NULL, **eigvec2 = NULL;
1066     matrix         topbox;
1067     real           xid, totmass, *sqrtm, *w_rls, t, lambda;
1068     int            natoms, step;
1069     char          *grpname;
1070     const char    *indexfile;
1071     char           title[STRLEN];
1072     int            i, j, d;
1073     int            nout, *iout, noutvec, *outvec, nfit;
1074     atom_id       *index, *ifit;
1075     const char    *VecFile, *Vec2File, *topfile;
1076     const char    *EigFile, *Eig2File;
1077     const char    *CompFile, *RmsfFile, *ProjOnVecFile;
1078     const char    *TwoDPlotFile, *ThreeDPlotFile;
1079     const char    *FilterFile, *ExtremeFile;
1080     const char    *OverlapFile, *InpMatFile;
1081     gmx_bool       bFit1, bFit2, bM, bIndex, bTPS, bTop, bVec2, bProj;
1082     gmx_bool       bFirstToLast, bFirstLastSet, bTraj, bCompare, bPDB3D;
1083     real          *eigval1 = NULL, *eigval2 = NULL;
1084     int            neig1, neig2;
1085     double       **xvgdata;
1086     output_env_t   oenv;
1087     gmx_rmpbc_t    gpbc;
1088
1089     t_filenm       fnm[] = {
1090         { efTRN, "-v",    "eigenvec",    ffREAD  },
1091         { efTRN, "-v2",   "eigenvec2",   ffOPTRD },
1092         { efTRX, "-f",    NULL,          ffOPTRD },
1093         { efTPS, NULL,    NULL,          ffOPTRD },
1094         { efNDX, NULL,    NULL,          ffOPTRD },
1095         { efXVG, "-eig", "eigenval",     ffOPTRD },
1096         { efXVG, "-eig2", "eigenval2",   ffOPTRD },
1097         { efXVG, "-comp", "eigcomp",     ffOPTWR },
1098         { efXVG, "-rmsf", "eigrmsf",     ffOPTWR },
1099         { efXVG, "-proj", "proj",        ffOPTWR },
1100         { efXVG, "-2d",   "2dproj",      ffOPTWR },
1101         { efSTO, "-3d",   "3dproj.pdb",  ffOPTWR },
1102         { efTRX, "-filt", "filtered",    ffOPTWR },
1103         { efTRX, "-extr", "extreme.pdb", ffOPTWR },
1104         { efXVG, "-over", "overlap",     ffOPTWR },
1105         { efXPM, "-inpr", "inprod",      ffOPTWR }
1106     };
1107 #define NFILE asize(fnm)
1108
1109     if (!parse_common_args(&argc, argv,
1110                            PCA_CAN_TIME | PCA_TIME_UNIT | PCA_CAN_VIEW | PCA_BE_NICE,
1111                            NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
1112     {
1113         return 0;
1114     }
1115
1116     indexfile = ftp2fn_null(efNDX, NFILE, fnm);
1117
1118     VecFile         = opt2fn("-v", NFILE, fnm);
1119     Vec2File        = opt2fn_null("-v2", NFILE, fnm);
1120     topfile         = ftp2fn(efTPS, NFILE, fnm);
1121     EigFile         = opt2fn_null("-eig", NFILE, fnm);
1122     Eig2File        = opt2fn_null("-eig2", NFILE, fnm);
1123     CompFile        = opt2fn_null("-comp", NFILE, fnm);
1124     RmsfFile        = opt2fn_null("-rmsf", NFILE, fnm);
1125     ProjOnVecFile   = opt2fn_null("-proj", NFILE, fnm);
1126     TwoDPlotFile    = opt2fn_null("-2d", NFILE, fnm);
1127     ThreeDPlotFile  = opt2fn_null("-3d", NFILE, fnm);
1128     FilterFile      = opt2fn_null("-filt", NFILE, fnm);
1129     ExtremeFile     = opt2fn_null("-extr", NFILE, fnm);
1130     OverlapFile     = opt2fn_null("-over", NFILE, fnm);
1131     InpMatFile      = ftp2fn_null(efXPM, NFILE, fnm);
1132
1133     bTop   = fn2bTPX(topfile);
1134     bProj  = ProjOnVecFile || TwoDPlotFile || ThreeDPlotFile
1135         || FilterFile || ExtremeFile;
1136     bFirstLastSet  =
1137         opt2parg_bSet("-first", NPA, pa) && opt2parg_bSet("-last", NPA, pa);
1138     bFirstToLast = CompFile || RmsfFile || ProjOnVecFile || FilterFile ||
1139         OverlapFile || ((ExtremeFile || InpMatFile) && bFirstLastSet);
1140     bVec2  = Vec2File || OverlapFile || InpMatFile;
1141     bM     = RmsfFile || bProj;
1142     bTraj  = ProjOnVecFile || FilterFile || (ExtremeFile && (max == 0))
1143         || TwoDPlotFile || ThreeDPlotFile;
1144     bIndex = bM || bProj;
1145     bTPS   = ftp2bSet(efTPS, NFILE, fnm) || bM || bTraj ||
1146         FilterFile  || (bIndex && indexfile);
1147     bCompare = Vec2File || Eig2File;
1148     bPDB3D   = fn2ftp(ThreeDPlotFile) == efPDB;
1149
1150     read_eigenvectors(VecFile, &natoms, &bFit1,
1151                       &xref1, &bDMR1, &xav1, &bDMA1,
1152                       &nvec1, &eignr1, &eigvec1, &eigval1);
1153     neig1 = DIM*natoms;
1154
1155     /* Overwrite eigenvalues from separate files if the user provides them */
1156     if (EigFile != NULL)
1157     {
1158         int neig_tmp = read_xvg(EigFile, &xvgdata, &i);
1159         if (neig_tmp != neig1)
1160         {
1161             fprintf(stderr, "Warning: number of eigenvalues in xvg file (%d) does not mtch trr file (%d)\n", neig1, natoms);
1162         }
1163         neig1 = neig_tmp;
1164         srenew(eigval1, neig1);
1165         for (j = 0; j < neig1; j++)
1166         {
1167             real tmp = eigval1[j];
1168             eigval1[j] = xvgdata[1][j];
1169             if (debug && (eigval1[j] != tmp))
1170             {
1171                 fprintf(debug, "Replacing eigenvalue %d. From trr: %10g, from xvg: %10g\n",
1172                         j, tmp, eigval1[j]);
1173             }
1174         }
1175         for (j = 0; j < i; j++)
1176         {
1177             sfree(xvgdata[j]);
1178         }
1179         sfree(xvgdata);
1180         fprintf(stderr, "Read %d eigenvalues from %s\n", neig1, EigFile);
1181     }
1182
1183     if (bEntropy)
1184     {
1185         if (bDMA1)
1186         {
1187             gmx_fatal(FARGS, "Can not calculate entropies from mass-weighted eigenvalues, redo the analysis without mass-weighting");
1188         }
1189         calc_entropy_qh(stdout, neig1, eigval1, temp, nskip);
1190         calc_entropy_schlitter(stdout, neig1, nskip, eigval1, temp);
1191     }
1192
1193     if (bVec2)
1194     {
1195         if (!Vec2File)
1196         {
1197             gmx_fatal(FARGS, "Need a second eigenvector file to do this analysis.");
1198         }
1199         read_eigenvectors(Vec2File, &neig2, &bFit2,
1200                           &xref2, &bDMR2, &xav2, &bDMA2, &nvec2, &eignr2, &eigvec2, &eigval2);
1201
1202         neig2 = DIM*neig2;
1203         if (neig2 != neig1)
1204         {
1205             gmx_fatal(FARGS, "Dimensions in the eigenvector files don't match");
1206         }
1207     }
1208
1209     if (Eig2File != NULL)
1210     {
1211         neig2 = read_xvg(Eig2File, &xvgdata, &i);
1212         srenew(eigval2, neig2);
1213         for (j = 0; j < neig2; j++)
1214         {
1215             eigval2[j] = xvgdata[1][j];
1216         }
1217         for (j = 0; j < i; j++)
1218         {
1219             sfree(xvgdata[j]);
1220         }
1221         sfree(xvgdata);
1222         fprintf(stderr, "Read %d eigenvalues from %s\n", neig2, Eig2File);
1223     }
1224
1225
1226     if ((!bFit1 || xref1) && !bDMR1 && !bDMA1)
1227     {
1228         bM = FALSE;
1229     }
1230     if ((xref1 == NULL) && (bM || bTraj))
1231     {
1232         bTPS = TRUE;
1233     }
1234
1235     xtop  = NULL;
1236     nfit  = 0;
1237     ifit  = NULL;
1238     w_rls = NULL;
1239
1240     if (!bTPS)
1241     {
1242         bTop = FALSE;
1243     }
1244     else
1245     {
1246         bTop = read_tps_conf(ftp2fn(efTPS, NFILE, fnm),
1247                              title, &top, &ePBC, &xtop, NULL, topbox, bM);
1248         atoms = &top.atoms;
1249         gpbc  = gmx_rmpbc_init(&top.idef, ePBC, atoms->nr);
1250         gmx_rmpbc(gpbc, atoms->nr, topbox, xtop);
1251         /* Fitting is only required for the projection */
1252         if (bProj && bFit1)
1253         {
1254             if (xref1 == NULL)
1255             {
1256                 printf("\nNote: the structure in %s should be the same\n"
1257                        "      as the one used for the fit in g_covar\n", topfile);
1258             }
1259             printf("\nSelect the index group that was used for the least squares fit in g_covar\n");
1260             get_index(atoms, indexfile, 1, &nfit, &ifit, &grpname);
1261
1262             snew(w_rls, atoms->nr);
1263             for (i = 0; (i < nfit); i++)
1264             {
1265                 if (bDMR1)
1266                 {
1267                     w_rls[ifit[i]] = atoms->atom[ifit[i]].m;
1268                 }
1269                 else
1270                 {
1271                     w_rls[ifit[i]] = 1.0;
1272                 }
1273             }
1274
1275             snew(xrefp, atoms->nr);
1276             if (xref1 != NULL)
1277             {
1278                 /* Safety check between selected fit-group and reference structure read from the eigenvector file */
1279                 if (natoms != nfit)
1280                 {
1281                     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);
1282                 }
1283                 for (i = 0; (i < nfit); i++)
1284                 {
1285                     copy_rvec(xref1[i], xrefp[ifit[i]]);
1286                 }
1287             }
1288             else
1289             {
1290                 /* The top coordinates are the fitting reference */
1291                 for (i = 0; (i < nfit); i++)
1292                 {
1293                     copy_rvec(xtop[ifit[i]], xrefp[ifit[i]]);
1294                 }
1295                 reset_x(nfit, ifit, atoms->nr, NULL, xrefp, w_rls);
1296             }
1297         }
1298         gmx_rmpbc_done(gpbc);
1299     }
1300
1301     if (bIndex)
1302     {
1303         printf("\nSelect an index group of %d elements that corresponds to the eigenvectors\n", natoms);
1304         get_index(atoms, indexfile, 1, &i, &index, &grpname);
1305         if (i != natoms)
1306         {
1307             gmx_fatal(FARGS, "you selected a group with %d elements instead of %d", i, natoms);
1308         }
1309         printf("\n");
1310     }
1311
1312     snew(sqrtm, natoms);
1313     if (bM && bDMA1)
1314     {
1315         proj_unit = "u\\S1/2\\Nnm";
1316         for (i = 0; (i < natoms); i++)
1317         {
1318             sqrtm[i] = sqrt(atoms->atom[index[i]].m);
1319         }
1320     }
1321     else
1322     {
1323         proj_unit = "nm";
1324         for (i = 0; (i < natoms); i++)
1325         {
1326             sqrtm[i] = 1.0;
1327         }
1328     }
1329
1330     if (bVec2)
1331     {
1332         t       = 0;
1333         totmass = 0;
1334         for (i = 0; (i < natoms); i++)
1335         {
1336             for (d = 0; (d < DIM); d++)
1337             {
1338                 t       += sqr((xav1[i][d]-xav2[i][d])*sqrtm[i]);
1339                 totmass += sqr(sqrtm[i]);
1340             }
1341         }
1342         fprintf(stdout, "RMSD (without fit) between the two average structures:"
1343                 " %.3f (nm)\n\n", sqrt(t/totmass));
1344     }
1345
1346     if (last == -1)
1347     {
1348         last = natoms*DIM;
1349     }
1350     if (first > -1)
1351     {
1352         if (bFirstToLast)
1353         {
1354             /* make an index from first to last */
1355             nout = last-first+1;
1356             snew(iout, nout);
1357             for (i = 0; i < nout; i++)
1358             {
1359                 iout[i] = first-1+i;
1360             }
1361         }
1362         else if (ThreeDPlotFile)
1363         {
1364             /* make an index of first+(0,1,2) and last */
1365             nout = bPDB3D ? 4 : 3;
1366             nout = min(last-first+1, nout);
1367             snew(iout, nout);
1368             iout[0] = first-1;
1369             iout[1] = first;
1370             if (nout > 3)
1371             {
1372                 iout[2] = first+1;
1373             }
1374             iout[nout-1] = last-1;
1375         }
1376         else
1377         {
1378             /* make an index of first and last */
1379             nout = 2;
1380             snew(iout, nout);
1381             iout[0] = first-1;
1382             iout[1] = last-1;
1383         }
1384     }
1385     else
1386     {
1387         printf("Select eigenvectors for output, end your selection with 0\n");
1388         nout = -1;
1389         iout = NULL;
1390
1391         do
1392         {
1393             nout++;
1394             srenew(iout, nout+1);
1395             if (1 != scanf("%d", &iout[nout]))
1396             {
1397                 gmx_fatal(FARGS, "Error reading user input");
1398             }
1399             iout[nout]--;
1400         }
1401         while (iout[nout] >= 0);
1402
1403         printf("\n");
1404     }
1405     /* make an index of the eigenvectors which are present */
1406     snew(outvec, nout);
1407     noutvec = 0;
1408     for (i = 0; i < nout; i++)
1409     {
1410         j = 0;
1411         while ((j < nvec1) && (eignr1[j] != iout[i]))
1412         {
1413             j++;
1414         }
1415         if ((j < nvec1) && (eignr1[j] == iout[i]))
1416         {
1417             outvec[noutvec] = j;
1418             noutvec++;
1419         }
1420     }
1421     fprintf(stderr, "%d eigenvectors selected for output", noutvec);
1422     if (noutvec <= 100)
1423     {
1424         fprintf(stderr, ":");
1425         for (j = 0; j < noutvec; j++)
1426         {
1427             fprintf(stderr, " %d", eignr1[outvec[j]]+1);
1428         }
1429     }
1430     fprintf(stderr, "\n");
1431
1432     if (CompFile)
1433     {
1434         components(CompFile, natoms, eignr1, eigvec1, noutvec, outvec, oenv);
1435     }
1436
1437     if (RmsfFile)
1438     {
1439         rmsf(RmsfFile, natoms, sqrtm, eignr1, eigvec1, noutvec, outvec, eigval1,
1440              neig1, oenv);
1441     }
1442
1443     if (bProj)
1444     {
1445         project(bTraj ? opt2fn("-f", NFILE, fnm) : NULL,
1446                 bTop ? &top : NULL, ePBC, topbox,
1447                 ProjOnVecFile, TwoDPlotFile, ThreeDPlotFile, FilterFile, skip,
1448                 ExtremeFile, bFirstLastSet, max, nextr, atoms, natoms, index,
1449                 bFit1, xrefp, nfit, ifit, w_rls,
1450                 sqrtm, xav1, eignr1, eigvec1, noutvec, outvec, bSplit,
1451                 oenv);
1452     }
1453
1454     if (OverlapFile)
1455     {
1456         overlap(OverlapFile, natoms,
1457                 eigvec1, nvec2, eignr2, eigvec2, noutvec, outvec, oenv);
1458     }
1459
1460     if (InpMatFile)
1461     {
1462         inprod_matrix(InpMatFile, natoms,
1463                       nvec1, eignr1, eigvec1, nvec2, eignr2, eigvec2,
1464                       bFirstLastSet, noutvec, outvec);
1465     }
1466
1467     if (bCompare)
1468     {
1469         compare(natoms, nvec1, eigvec1, nvec2, eigvec2, eigval1, neig1, eigval2, neig2);
1470     }
1471
1472
1473     if (!CompFile && !bProj && !OverlapFile && !InpMatFile &&
1474         !bCompare && !bEntropy)
1475     {
1476         fprintf(stderr, "\nIf you want some output,"
1477                 " set one (or two or ...) of the output file options\n");
1478     }
1479
1480
1481     view_all(oenv, NFILE, fnm);
1482
1483     return 0;
1484 }