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