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