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