src/gromacs/gmxana/gmx_anaeig.c

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